Abstract
The ongoing transition to renewable energy supply comes with a restructuring of power grids, changing their effective interaction topologies, more and more strongly decentralizing them and substantially modifying their input, output, and response characteristics. All of these changes imply that power grids become increasingly affected by collective, nonlinear dynamic phenomena, structurally and dynamically more distributed and less predictable in space and time, more heterogeneous in its building blocks, and as a consequence less centrally controllable. Here cornerstone aspects of data-driven and mathematical modeling of collective dynamical phenomena emerging in real and model power grid networks by combining theories from nonlinear dynamics, stochastic processes and statistical physics, anomalous statistics, optimization, and graph theory are reviewed. The mathematical background required for adequate modeling and analysis approaches is introduced, an overview of power system models is given, and a range of collective dynamical phenomena are focused on, including synchronization and phase locking, flow (re)routing, Braess’s paradox, geometric frustration, and spreading and localization of perturbations and cascading failures, as well as the nonequilibrium dynamics of power grids, where fluctuations play a pivotal role.
27 More- Received 30 April 2021
DOI:https://doi.org/10.1103/RevModPhys.94.015005
© 2022 American Physical Society
Physics Subject Headings (PhySH)
- Physical Systems
Authors & Affiliations
- Forschungszentrum Jülich, Institute for Energy and Climate Research (IEK-STE), 52428 Jülich, Germany and Institute for Theoretical Physics, University of Cologne, Köln 50937, Germany
- Potsdam Institute for Climate Impact Research, Potsdam 14473, Germany
- Potsdam Institute for Climate Impact Research, Potsdam 14473, Germany, Humboldt University of Berlin, Berlin 12489, Germany, and Lobachevsky University of Nizhny Novgorod, Nizhnij Novgorod 603950, Russia
- Department of Physics and Earth Sciences and Department of Computer Science, Jacobs University Bremen, Campus Ring 1, 28759 Bremen, Germany and Division of Advanced Materials Science, POSTECH, San 31, Hyoja-dong, Nam-gu, Pohang 790-784, Republic of Korea
- Department of Physics and Earth Sciences, Jacobs University Bremen, Campus Ring 1, 28759 Bremen, Germany
- Chair for Network Dynamics, Center for Advancing Electronics (cfaed) and Institute for Theoretical Physics, Technical University of Dresden, 01062 Dresden, Germany and Lakeside Labs, 9020 Klagenfurt am Wörthersee, Austria
- *d.witthaut@fz-juelich.de
- †hellmann@pik-potsdam.de
- ‡kurths@pik-potsdam.de
- §s.kettemann@jacobs-university.de
- ¶h.ortmanns@jacobs-university.de
- ∥marc.timme@tu-dresden.de
Article Text
I. Introduction
The mitigation of climate change is one of the greatest challenges of mankind. Currently, 65% of all greenhouse gas emissions are caused by the carbon dioxide ( CO2) emissions from fossil fuel combustion and industrial processes (IPCC, 2014). Thus, a fundamental transformation of our energy system is inevitable to meet the goals of the Paris agreement to limit global warming (Rogelj et al., 2015). Rapid actions are needed as greenhouse gas emissions continue to increase and the remaining carbon budgets shrink at an alarming rate (Figueres et al., 2017; Rockström et al., 2017).
Power grids are at the heart of this transformation. Renewable energy sources have shown remarkable development in recent decades (Edenhofer, Pichs-Madruga, and Sokona, 2011; Wiser et al., 2016; Creutzig et al., 2017), but they fundamentally change the operation of the grids they are connected to. The decarbonization of heating, transport, and other sectors introduces new consumers to the electric power grid. Gas and oil heating is replaced by heat pumps, and district heating plants use resistive heaters in times of high renewable power generation (Gröger, Gasteiger, and Suchsland, 2015). Batteries of electric vehicles have to be charged, which may cause congestion in the grid (Carvalho et al., 2015; Bloess, Schill, and Zerrahn, 2018) but also offer a chance of improving stability (Liu et al., 2013; Gajduk et al., 2014).
Power grids are central to our economy and daily life. It is a “uniquely critical infrastructure,” as a variety of other infrastructures and sectors are dependent on the secure supply of electric power (Van der Vleuten and Lagendijk, 2010). Communication, information processing, public transport, and even cooling infrastructure for food supply cannot work without electric power. The energy transition further increases this dependence, as other sectors are being electrified. Ensuring a stable operation of the power system is thus a necessity and a grand challenge.
Transforming the electric power system to meet the 21st century challenges constitutes a transdisciplinary challenge (Brummitt et al., 2013). Modeling, simulation, and analysis of the physical grid have always been an integral part of power engineering. By now, thousands of simulations are run every hour to assess the state of the grid and potential contingencies. But the grid does not run in isolation: it is affected by the weather, by automatic controls and information systems, by energy markets, and finally also by human operators and consumers. Understanding such interdependent systems is challenging, particularly when they leave the normal state of operation in the case of contingencies. Complexity science can effectively complement detailed simulation models. It can elucidate fundamental mechanisms and interactions, provide explicative models, and identify new and unexpected risks. At the same time large amounts of data become available from new measurement devices and smart grid infrastructures. This fosters an empirical and statistical view on power grid operation and stability. In summary, power systems are becoming increasingly complex and methods from statistical physics, nonlinear dynamics, and network science can play an important part in addressing current challenges.
We now name three central challenges for power system operation arising during the energy transition:
(i) The variability of renewable energy sources is perhaps the most well-known one. Wind and solar power generation is determined by the weather and fluctuates on all temporal and spatial scales. Synoptic and seasonal variability are the most obvious features (Heide et al., 2010; Bloomfield et al., 2018; Staffell and Pfenninger, 2018), but power generation can also fluctuate within seconds due to atmospheric turbulence and the dynamics of clouds (Milan, Wächter, and Peinke, 2013; Anvari, Wächter, and Peinke, 2016; Anvari et al., 2016; Zhang et al., 2019). Fingerprints of these fluctuations can be observed in power grid frequency data (Haehne et al., 2018); cf. also Gorjão, Jumar et al. (2020) for a comprehensive data source.
Longer timescales have attracted increasing interest recently, ranging from interannual (Collins et al., 2018) to decadal variability (Wohland et al., 2019) to the impact of climate change (Wohland et al., 2017; Schlott et al., 2018; Weber et al., 2018). A thorough statistical understanding of all modes of variability and their consequences is integral to system operation.
(ii) Furthermore, wind and solar power are typically generated at locations with favorable natural resources. These locations are often far from consumers, so a long-distance transmission of electric power becomes necessary (Pesch, Allelein, and Hake, 2014). This increases grid loads and makes it vulnerable to large-scale blackouts via cascading failures. Robustness and vulnerability are central topics in network science. Physicists have developed and analyzed a variety of models for cascading failures to understand their propagation, their statistical features, and methods to mitigate them; see Albert, Jeong, and Barabási (2000), Carreras et al. (2002), Motter and Lai (2002), Dobson et al. (2007), Witthaut et al. (2016), Yang, Nishikawa, and Motter (2017), Nesti, Zocca, and Zwart (2018), and Schäfer, Witthaut et al. (2018). More recently the optimal design of networks has attracted strong interest in the scientific community; see Nishikawa and Motter (2006), Katifori, Szöllősi, and Magnasco (2010), and Kaiser, Ronellenfitsch, and Witthaut (2020).
(iii) Finally, renewable power sources are different from conventional power plants operating synchronous machines. They are typically connected to the grid via power electronic inverters with different characteristics and dynamics (Carrasco et al., 2006). Inverters have some important disadvantages, in particular, they have no intrinsic inertia to stabilize grid dynamics (Milano et al., 2018) but offer great flexibility in design and control. Hence, the dynamical stability of complex networked systems is a topic of rapidly growing importance; cf. Anvari, Hellmann, and Zhang (2020). Research must address the stability of existing systems as well as the design of future systems.
The goal of this review is twofold. It first provides a starting point for physicists interested in aspects of power system dynamics, operation, and robustness. To this end, the first part of the review is written as a tutorial. We provide a short review of the main tools from dynamical systems and network science in Sec. II. We then review basic principles of power grid operation and provide an overview of static and dynamic models for power grids in Sec. III. We focus on the mathematical description of the elements of power systems, avoiding most technical details covered in the engineering literature (Grainger and Stevenson, 1994; Kundur, 1994; Machowski, Bialek, and Bumby, 2008; Wood, Wollenberg, and Sheblé, 2014). Section IV then discusses the availability of power system datasets required for own studies and simulations. The tutorial part closes with a review of fundamental aspects of power grid stability in Sec. V.
The second part of the review addresses some recent results obtained at the interface of statistical physics, dynamical systems, network science, and power engineering. We first consider aspects of dynamical stability in Sec. VI. When do stable fixed points exist and how do they depend on local dynamics and network topologies? We then analyze the robustness and vulnerability of power grids in Sec. VII. How does a grid react to the failure of single elements and how do large-scale blackouts emerge? We investigate dynamical perturbations or fluctuations and the transient dynamics of grids in Sec. VIII. We conclude with an outlook on current and future research topics in Sec. IX.
II. Fundamentals of Graphs and Networks
A broad variety of collective dynamical phenomena in future-compliant power systems crucially depend on how consumers, producers, and storage infrastructures are interconnected. The topology of the resulting interactions typically exhibits a variety of complex structural features beyond regular lattices or random structures. We introduce here basic concepts from graph theory used to quantitatively describe such topologies and remark by example where they emerge in the analysis and modeling of power systems. A more comprehensive introduction to graph theory and its applications can be found in a variety of textbooks, such as Bollobás (1998).
A graph G=(V,E) is given by a set V={1,…,N} of vertices together with a set E of edges that is a symmetric subset of V×V; i.e., if i,j∈V then (i,j)∈E if and only if (j,i)∈E. Because of this symmetry, such a graph is also called an undirected graph. In a directed graph, the edges are ordered 2-tuples, denoted (i,j)∈E, thus have a direction from j to i and the edge set is not necessarily symmetric.
If (i,j)∈E, we call i and j adjacent or neighbors and write i∼j, where for directed graphs i might be adjacent to j but not vice versa. The degree of a vertex is the number of its neighbors; for directed graphs the in-degree is the number of edges pointing to the vertex, and the out-degree is the number of edges pointing from the vertex. If the degree of every vertex is the same, ki≡k for all i, the graph is called k-regular or simply regular.
A walk in a graph is a sequence of vertices and edges
where each successive pair of vertices vj,vj+1 is connected by an edge ej=(vj,vj+1). If no vertex repeats in such a sequence, it is called a path. A cycle is a walk with at least one edge, with first and last vertices identical and with no other vertices repeating. A graph is connected if there is a path from any of its vertices to any other. A directed graph is strongly connected if there is a directed path from any of its vertices to any other; it is weakly connected if for any two of its vertices v and w there is a directed path from v to w or from w to v. A graph with no cycles is called a forest; a connected forest is a tree.
Paths can be used to unambiguously define a metric on a graph. The unweighted length of a path p is simply the number k of edges in the sense of Eq. (1). The geodesic distance between two vertices vn and vm is defined as the length of the shortest path connecting the two vertices. Note that the shortest path is not necessarily unique; if no such path exists, the distance is infinite. In many applications, such a simple definition of distance is not sufficient, as edges may exhibit heterogeneous features. If a distance or weight d(e) is assigned to every edge e∈E, the weighted length of a path is given by
As before the length defines a weighted geodesic distance via the shortest path, provided that all d(e) are non-negative.
The modeling of power grids and other real-world networks includes but goes far beyond the basic concepts of graph theory. In particular, we may assign quantities such as voltages or power injections to the nodes, currents, and power flows to the edges and study dynamical systems on networks. We will thus need an algebraic description of networks, which links properties of nodes and edges to the previously defined basic properties of graphs. For a more detailed overview of algebraic graph theory and its applications to power grids, see Dörfler, Simpson-Porco, and Bullo (2018).
A power grid has N vertices, often referred to as buses in power engineering, and L edges or branches. Properties of vertices such as voltages are encoded in vectors in RN or CN. The absolute value of each complex number encodes the voltage magnitude and the arguments encode its phase (relative to a reference). Similarly, properties of edges such as flows are encoded in vectors in RL or CL. The topology of such a network or graph is encoded in the adjacency matrix A∈RN×N with the components
Furthermore, it is useful to define the node-edge incidence matrix E∈RN×L with the components (Newman, 2010)
We note that power grids abstract to undirected graphs: a transmission line can transmit power in either direction. For each edge we fix an orientation to encode the direction of power flow: positive in one direction and negative in the other one. The orientation, the start, and the end of an edge is arbitrary but must be kept fixed.
The node-edge incidence matrix is useful to keep track of network flows. Suppose that the flows over all edges are denoted as a vector F∈RL. The total inflow or outflow at all nodes is then given by EF∈RN. The kernel of E corresponds to cycle flows: For any vector F≠0 in the kernel, the inflow and outflow vanish at every node such that the flow must be cyclic. We can fix a basis for the kernel by choosing L−N+1 independent fundamental cycles and encoding this basis in the cycle-edge incidence matrix C∈RL×(L−N+1) with components
such that EC=0. In a plane graph, a graph drawn on a plane without edge crossings, one typically chooses the facets of the graph as the fundamental cycles. After fixing a basis, any cycle flow (any vector in the kernel of E) can be written as F=Cf, with coefficients f1,…,fL−N+1.
Not all edges in a grid are equally strong; for instance, the impedance typically scales with the length of a line such that longer transmission lines will have lower conductances and susceptances. We assign weights bij to all edges (i,j). Either these weights can be included in a diagonal matrix
or we can define a weighted adjacency matrix
A particularly important matrix to characterize the stability of a network dynamical system is the Laplacian matrix Λ∈RN×N with components (Newman, 2010)
Using the node-edge incidence matrix, the Laplacian is written as Λ=EBdE⊤, where the superscript ⊤ denotes the transpose of a matrix. If all weights are positive, Λ is symmetric and positive semidefinite; i.e., all eigenvalues are real and non-negative and can be ordered as
The second eigenvalue λ2 measures the algebraic connectivity of a graph (Fiedler, 1973), and the associated eigenvector v2 is called the Fiedler vector. More generally, the number of zero eigenvalues of Λ equals the number of connected components of the graph (Newman, 2010).
III. Power System Models
In this section we introduce the basic mathematical models for power system operation and stability. We start with the fundamental relations for currents and power flows in ac power grids leading to static power grid models, then introduce dynamical models for synchronous machines and power electronic inverters.
A. Static models
1. Fundamentals of ac power flow
In an ac power grid the voltages V and currents I oscillate approximately sinusoidally with time as follows:
where we have defined the complex-valued amplitudes as
Voltage and current are typically out of phase due to capacities or inductances, which is reflected by the phase factors ϑV and ϑI. As the phases of voltage and current play essential roles in describing grid operation, we mainly employ a complex notation in this review. In the engineering literature, the complex quantities V_ and I_ are referred to as phasors and denoted by an underline. We adopt this notation in the following but denote the imaginary unit by the symbol i, as is common in physics. The electric power is oscillating too, but only the time-averaged value
the real power or active power, can do work. The imaginary part Q=I(V_I_∗) describes the power temporarily stored in capacitances and inductances and is referred to as the reactive power. Furthermore, one defines the apparent power S=P+iQ=V_I_∗. These relations show that the relative phases of voltages and currents are essential for the power flow in a grid. We note that deviations from perfect sinusoidal signals do occur in practice: for example, higher harmonics in power electronic sources (Liang and Andalib-Bin-Karim, 2018). Nevertheless, this idealization is extremely helpful for the modeling and analysis of power systems.
The basic relation between voltage at the nodes k,n∈{1,…,N} of a power grid and the currents flowing between the nodes is given by Ohm’s law:
where zkn=rkn+ixkn is the impedance and ykn=1/zkn is the admittance of the transmission line between the nodes k and n. The admittance is divided into its real and imaginary parts ykn=gkn+ibkn, where gkn is the conductance and bkn is the susceptance. For multicircuit transmission lines, we take ykn to be the sum of the admittances of the single circuits and ykn=0 if no transmission line exists. The total current injected into the node k is then given by I_k=∑nI_kn=∑nykn(V_k−V_n). For actual calculations it is convenient to introduce the nodal admittance matrix Y∈CN×N with the entries
Ohm’s law can then be rewritten in a vectorial form, which yields the network equations
where I_=(I_1,I_2,…,I_N) represents the complex currents injected into the N nodes and V_=(V_1,V_2,…,V_N) denotes the complex voltages. The apparent power injected to a node k is then given by
Real-world power transmission and distribution grids are mostly constructed as three-phase systems. Hence, there are three conductors with voltages VA(t)=Vpeakcos(ωt+ϑV), VB(t)=Vpeakcos(ωt+ϑV+2π/3), and VC(t)=Vpeakcos(ωt+ϑV+4π/3) to ground. The voltage between two conductors oscillates sinusoidally with an amplitude of √3Vpeak. If all three conductors have the same voltage magnitude and a fixed phase difference, it is sufficient to keep one value of V_ and I_. Basic relations such as Eq. (15) remain unchanged, but the transmitted power is 3 times as large as in a single-phase system. The power balance at each node hence reads
where θkn=θk−θn is the phase difference of complex voltages V_k and V_n. Real and imaginary parts yield the active and the reactive power balance conditions
Detailed discussions of the network equations and the resulting power flow equations were provided in Chaps. 7–9 of Grainger and Stevenson (1994) and Chap. 6 of Wood, Wollenberg, and Sheblé (2014).
2. Transformers and the per unit system
Real-world transmission systems are more diverse and include different transmission elements. In particular, we consider the following three types of devices. First, transformers are essential to link different voltage levels, from the transmission grid at highest voltages (380 kV in Europe, up to 765 kV in North America) down to the house connecting lines at low voltages (380 V in Europe). Second, special transformers may shift the phase of voltage and current on the two terminal ends (Verboomen et al., 2005). Third, the model for an ordinary transmission line needs to be extended to take into account the charging capacity typically present in real-world lines.
One can include all these devices in the common network model using a system of rescaled variables, the per unit (pu) system, and the unified model for all transmission elements depicted in Fig. 1(b). The equivalent circuit of this transmission element includes an ideal transformer with the tap ratio t and the phase shift θshift in addition to a π-equivalent line model. The line model includes the series admittance y and a charging susceptance bc, which is attributed equally to the two end points of the line for simplicity.
To see how the network equations have to be modified in each case, we first consider a single transmission element as shown in Fig. 1(b). The voltages and currents on the two terminals of the ideal transformer are related by a factor of N12=t12eiθshift12. The current flowing through the series admittance is given by y12(V_2−V_1/N12) according to Ohm’s law. Kirchhoff’s current law for the junctions marked by thick black dots in the figure reads
It is convenient to introduce scaled units, which are referred to as the “per unit” or “pu” system in power engineering. We fix a reference value Sbase for the power in the entire grid and a reference value Vbase,k for the voltage separately for every voltage level k. The reference values for the currents and admittances are then given by Ibase,k=Sbase/(3Vbase,k) and Ybase,k=Ibase,k/Vbase,k. Generally, one selects the nominal voltage of each voltage level as the reference value, while a typical value for the power would be Sbase=100 MVA. For the π-transmission line depicted in Fig. 1, this choice yields Vbase,1/Vbase,2=t12. Dividing Eq. (20) by Ibase,2=Ybase,2Vbase,2 then yields
where we use a tilde to denote the scaled quantities in the pu system, as with ˜I_2=I_2/Ibase,2. Equation (21) can be rescaled analogously if we adopt the convention that θshift12=−θshift21. We note that Ybase,2 and not Ybase,1 has to be used in the normalization of the impedance ˜y12=˜y21=y12/Ybase,2, as the line is on the 2 side of the transformer, not on the 1 side.
In a large power grid with many nodes and transmission lines, we then obtain
We note that an ideal transformer appears as a simple series admittance ˜y12 in the pu system. This simplification will be used in the modeling of generators, which are typically connected to the grid via a step-up transformer.
Rewriting Eq. (23) in a vectorial form, we recover the network equations in the form of Eq. (15), however, with a modified nodal admittance matrix ˜Y defined as
Finally, the power is rescaled to the pu system by dividing Eq. (17) by the reference value Sbase=3Vbase,kIbase,k, which yields ˜Sk=˜V_k˜I_∗k as desired.
The pu system thus simplifies the network equations, as we can treat all transmission elements alike. In addition, it is advantageous for numerical calculations, as all quantities are of order unity such that rounding errors are minimized. All details of the transmission elements are absorbed into the nodal admittance matrix ˜Y. In the absence of phase shifters, ˜Y=˜Y⊤ is symmetric. A detailed discussion of the pu systems and their benefits was provided in Chaps. 1 and 2 of Grainger and Stevenson (1994).
3. The ac load flow equations
Thus far we have established a set of nonlinear algebraic equations describing the steady state of ac power grids, linking nodal power injections Pn and Qn to nodal state variables θn and |Vn|. Two of these four quantities are fixed externally for every node, for instance, via power demand. The remaining unknown variables are obtained by numerical solution of the load flow equations; see Chap. 9 of Grainger and Stevenson (1994) and Chap. 6 of Wood, Wollenberg, and Sheblé (2014) for detailed introductions.
In practice, we distinguish three types of nodes. A PV bus is typically connected to a generator, which provides a fixed power output P at fixed voltage magnitudes |V|. A PQ bus represents a node with a given consumption, such that the net injected active power P and reactive power Q are given, while the voltages at these nodes are unknown. A special kind of bus is the slack bus, which captures the supply of the power necessary to have overall power balance in the system to ensure that there is a steady-state solution. It acts as an ideal voltage source (fixed V_) where the parameters P and Q remain unspecified to balance active and reactive power during the iteration toward the steady-state solution. This is necessary because power losses on the transmission lines are not known a priori before the solution is obtained.
Suppose that we are interested in a transmission grid of Nℓ loads, Ng generators, plus one generator node taken as slack. We then have 2Nℓ+Ng unknown state variables: the phases of all Ng+Nℓ nonslack nodes and the voltage magnitudes of Nℓ load nodes. These unknown variables are fixed by 2Nℓ+Ng nonlinear algebraic equations as follows:
where n labels all nonslack nodes, while k labels load nodes only. The superscript “sp” indicates that values on the right side are specified beforehand and the functions Pn(⋅) and Qk(⋅) are given by Eqs. (18) and (19).
A common and effective way to solve this equation system is the Newton-Raphson method [cf. Chap. 9 of Grainger and Stevenson (1994) and Chap. 6 of Wood, Wollenberg, and Sheblé (2014)], which iteratively updates the state vector from an initial guess at the solution of Eqs. (25). The solution for the state vector describes the power grid in steady-state operation. We note that a set of nonlinear equations may exhibit a single well-defined solution but may also have no or several solutions. The lack of a solution indicates an unstable situation that is further discussed in Sec. V.A, while multistability is addressed in Sec. VI.C.3. Recent advances and challenges in numerical methods for ac power flow computations were discussed by Trias (2012) and Mehta, Molzahn, and Turitsyn (2016).
4. Linearized power flow
For small loads and losses in a power grid, load flow calculations simplify considerably [see Chap. 6.18 of Wood, Wollenberg, and Sheblé (2014)] via linearized power flow equations that are based on three approximations. (i) For each transmission element, the resistance and the charging susceptance are small compared to the reactance and are thus neglected. Hence, the admittance is purely imaginary [ ynk=1/(ixnk)]. (ii) Variations of the voltage magnitude are typically small in transmission grids, such that we can fix them at the reference value of the respective voltage level. In the pu system we thus write V_n=eiθn, and the power balance at a node n reads [cf. Eqs. (18) and (20)]
Technically, all nodes must then be considered PV buses, such that no equations for the reactive power must be taken into account. Taking the real part of Eq. (26), we obtain the following balance equation for the real power:
(iii) Finally, small loads imply that the phase differences across a transmission line are small such that the sine function is approximated to first order. The load flow calculations then reduce to the following set of linear equations:
where Pshiftn=∑kBnkθshiftnk accounts for the effects of the phase shifting transformers and the matrix B∈RN×N is the nodal susceptance matrix expressed in per units but without taking into account potential phase shifts:
The simplified Eq. (28) is often referred to as the dc approximation, as it is mathematically equivalent to Kirchhoff’s circuit equation for dc electric circuits. Still, it describes the flow of real power in ac power grids. Linear equations can be solved much faster than the nonlinear load flow equations (25), which is advantageous when the flow must be calculated for many different scenarios of generation and load. Furthermore, the linear approximation avoids the problem of multiple or disappearing solutions or multiple optima in optimization problems. Limits of its applicability were discussed in detail in (Purchala et al. (2005), Van Hertem et al. (2006), and Stott, Jardim, and Alsac (2009).
A word of caution is in order, as the symbol B is used for several related but different quantities. Transmission lines and ordinary transformers have a positive reactance ˜x>0. Hence, the off-diagonal elements of the nodal susceptance matrix (29) are negative. In contrast, the imaginary parts of the nodal admittance matrix elements (14) are positive. Both quantities are denoted by the symbol B, and one must be careful to not confuse them.
5. Matrix formulation
The linearized power flow equations can be condensed in a highly practical compact matrix notation. Let N denote the number of nodes and L denote the number of lines in a grid as before. We define the vectors of power injection as P=(P1,…,PN)⊤∈RN, the vector of voltage phase angles as θ=(θ1,…,θN)⊤∈RN, and the vector of line flows as F=(F1,…,FL)⊤∈RL. As the flow is directed, we must assign an orientation to each line in the grid that is arbitrary but must be kept fixed. We assume that there are no phase-shifting elements and that power injections are balanced; i.e., the Pn sum to zero.
The directed real power flow Fℓ on a line ℓ from node m to node n is given by Fℓ=x−1ℓ(θm−θn). Using the diagonal matrix Bd defined in Eq. (6) with weights bℓ=x−1ℓ and the node-edge incidence matrix (4), the relation of flows and phase angles are given by
The real power balance at each node now reads
stating that the power flowing out of each node must equal the power injected at the node. When one combines Eqs. (30) and (31), the equation for the power injections in terms of the voltage angles is obtained as
When Eq. (32) is taken together with Eq. (30), we thus have a linear relation between line flows F and nodal power injections P.
The matrix B=EBdE⊤ with the elements given in Eq. (29) is a weighted Laplacian and thus has a single zero eigenvalue if the network is connected; cf. Sec. II. Thus, it is noninvertible and we solve Eq. (32) for θ via the Moore-Penrose pseudoinverse B∗ to obtain the line flows directly as a linear function of the nodal power injections
The matrix prefactor is called the power transfer distribution factor (PTDF) matrix (Wood, Wollenberg, and Sheblé, 2014), as it describes how power injections are distributed throughout the grid.
The zero eigenvalue of the Laplacian B represents a phase-shift invariance. Uniformly shifting all nodal phases by the same constant c, θn→θn+c, does not affect any power flows. Fixing the phase angle θk≡0 at one node k (the slack node) removes this degree of freedom. We then remove the node from the analysis by removing the kth row from the vector θ and the kth row and kth column from the matrix B. The resulting matrix is called a grounded Laplacian.
6. Generalized linear approximations
Several generalizations of the dc approximation have been developed, and we comment upon three of them.
(i) First, one can improve on the linear approximation of the sine function by rewriting the governing equations of the dc approximation in two parts. We have
where ψℓ denotes the sine of the phase difference along the line ℓ. In the ordinary dc approximation, one simply replaces the sine with its argument; in vectorial form one thus obtains ψ=E⊤θ. If we do not neglect the sine, we instead have
The basic idea of the generalized dc approximation (Dörfler and Bullo, 2013b; Simpson-Porco, 2018) is to first obtain ψ approximately and then solve for θ. The general solution to the underdetermined equation (34) reads
where C is the cycle-edge incidence matrix (5) and f is a vector of cycle flows. To find the actual value of the cycle flows, we would obtain the correct solution of the nonlinear power flow equations, an approach discussed in Sec. VI.C. Here we consider practical approximate solutions and assume the cycle flows f to be negligible under normal operating conditions. Equation (35) then yields
(ii) A second generalization is to approximately incorporate Ohmic losses in a linear fashion. Two iterative procedures that introduce losses but keep the assumption of fixed voltage magnitudes were described by Stott, Jardim, and Alsac (2009) and Simpson-Porco (2018).
(iii) A third class of linear models includes reactive power flows and voltage variations by appropriate linearization of the nonlinear ac load flow equations. Different approaches were discussed by Coffrin, Van Hentenryck, and Bent (2012), Zhang et al. (2013), and Yang et al. (2019).
7. Hybrid power grids
High voltage dc (HVDC) transmission lines can transmit bulk power over large distances with lower capital cost and lower losses than standard ac lines (Setreus and Bertling, 2008). They are connected to ac via power electronic converter stations, which offers a high degree of flexibility. The transmitted real power can be controlled by the grid operator. Modern converter stations can provide system services such as reactive power compensation. HVDC lines are included in a load flow study rather easily by representing the two converter stations with two additional PV buses. One of the converter nodes draws a power Pf<0 from the grid at node f, and the other converter then feeds a power
back into the grid at node t. Losses are typically small such that l0 is of the order of a few megawatts and l1 is of the order of a few percent.
B. Optimal power flow
A fundamental problem in power engineering and energy economics is to determine the cost-optimal dispatch of generators: Which generators should run to satisfy a given demand at minimum costs [see Chap. 13 of Grainger and Stevenson (1994) and Chap. 8 of Wood, Wollenberg, and Sheblé (2014)]? We address this problem here in terms of the dc approximation, taking into account the real power generation and demand but neglecting losses. Suppose that the total demand or load is given by Ploadtot and that there are Ng generators that can be used to satisfy this demand at minimum costs. The actual generation Pgenm of each generator m=1,…,Ng is bounded by the technical capacity such that we have the constraints
The power balance condition in the grid reads ∑Ngm=1Pgenm=Ploadtot. We assume for simplicity that the costs are proportional to the output Pgenm for each generator. The total variable costs can then be written as
where cm denotes the variable costs for generator m. If there are no further constraints, the solution with minimum costs is simple: Switch on all generators consecutively according to their variable costs cm such that
where c∗ is the variable cost of the last generator switched on that is identified with the market electricity price. The primed sum in Eq. (41) runs over all generators with cm<c∗. This is called the merit order principle.
In real-world applications there are many more necessary conditions and constraints. The real power flow in each transmission line must not exceed the line rating in order not to become overloaded and fail:
To satisfy this condition we must consider where power is generated and consumed. In particular, the demand or load Ploadi must be specified separately for every node i∈{1,…,N} of the grid. We can then express all line flows Fℓ in terms of the dispatch Pgenm as follows. We first introduce the following generator incidence matrix Egen, which indicates where each generator is connected to the grid:
The resulting net real power injection is thus given by
for every node n. Next, the line flows Fℓ can be expressed by power injections via Eq. (33), leading to an affine linear relation of flows Fℓ and optimization variables Pgenm. Finally, we arrive at the following important optimization problem, which is commonly referred to as dc or linear optimal power flow (OPF):
Mathematically, this optimization problem is a linear program that admits an efficient solution. Different formulations and their computational efficiency were discussed by Hörsch, Ronellenfitsch et al. (2018). Equation (45) represents the basic optimal power flow problem, and various extensions have been discussed for real-world applications leading to a vast body of literature (Huneault and Galiana, 1991; Momoh, Adapa, and El-Hawary, 1999; Momoh, El-Hawary, and Adapa, 1999; Frank, Steponavice, and Rebennack, 2012a, 2012b). The linear function (40) is a strong simplification of the real cost function. It can be replaced by a convex function without changing the overall complexity of the optimization problem. In contrast, a nonconvex cost function strongly alters the nature of the problem, making it hard to solve in general. Another approach is to consider different modes of operation of a generator. In the simplest case the generator is either on or off. If a generator is on, it typically has a minimum generation such that an additional constraint Pminm≤Pgenm arises. To take this switching behavior into account, an additional binary optimization variable is introduced, leading to a mixed-integer linear program.
Variable renewable energy sources may be included in two different ways. If renewable generation is taken to be fixed, the real power load at each node Ploadi is replaced by the residual load, i.e., the difference between the load and renewable generation. However, renewable generation will exceed the load at least temporarily in future energy systems. To allow for a curtailment of the power generation, wind farms and photovoltaic parks are included as generators with zero variable costs and Pmaxm is set to the assumed or forecasted power generation. For optimization, we may then choose any value Pgenm between 0 (complete curtailment) and Pmaxm, depending on what most benefits the grid. A systematic way of optimizing the dispatch under uncertainties in production and demand is sketched in Sec. VIII.E, where the system constraints are satisfied only with a specified probability.
A further natural extension is to abandon the assumptions of the dc approximation and start with the full nonlinear ac power flow equations. However, the resulting optimization problem is again nonconvex, which makes it much harder to solve. The development of algorithms for this problem is an active field of research; see Erseghe (2014) and Engelmann et al. (2017).
C. Dynamic models
While load flow calculations describe the steady state of a power grid, a large set of models of different complexity is available to analyze its dynamics and the stability of steady states. Which model to choose depends crucially on the phenomena and the timescales to be investigated (Sauer, Pai, and Chow, 2018). The bulk power generation and demand change on timescales of minutes to hours. Optimization models are routinely used to model how the unit commitment and power flows change from hour to hour. In the following, we focus on the next timescale and address aspects of synchronization, transient stability, and, to a lesser extent, voltage dynamics and load-frequency control. These phenomena typically take place on timescales of 10−1−101 s. Faster phenomena, including subtransient effects and electromagnetic propagation effects as well as the modeling of generator control and protection systems, are not covered in this review. Much more detailed and extensive introductions were given in standard textbooks (Kundur, 1994; Machowski, Bialek, and Bumby, 2008; Sauer, Pai, and Chow, 2018). Overviews of different modeling aspects were given by Gajduk, Todorovski, and Kocarev (2014) and Nishikawa and Motter (2015).
1. The swing equation
The classic swing equation describes the dynamics of the mechanical rotation of a synchronous machine; see Chap. 5.1 of Machowski, Bialek, and Bumby (2008) and Nishikawa and Motter (2015). The basic dynamical variable is the mechanical phase angle, which is identical to the voltage phase angle (hence the term “synchronous machine”). The phase angle δ is commonly measured relative to a frame of reference rotating at the reference frequency of the grid ωR. The dynamics is then simply given by Newton’s equation as
where ω=˙δ is the angular frequency, Δω=ω−ωR is its deviation from the reference, J is the moment of inertia of the machine, Tmech is the mechanical torque driving the rotating machine, and Tel is the electromagnetic torque due to the power transferred to the grid. The machine experiences mechanical friction with damping coefficient Dmech and an effective damping due to damper windings characterized by coefficient Del. In addition, one sometimes includes the effects of a frequency controller with gain constant κ, which is discussed in Sec. III.D. Aggregating the damping factors D=Dmech+Del+κ and noting that ˙δ=ω, we obtain
We have thus related the dynamics to the net mechanical input power Pmech and the electric power Pel acting on the machine. The swing equation is used to analyze the stability and small-amplitude dynamics around a steady state; hence, we can approximate ωR/ω≈1.
As for static flows, calculations are often carried out in dimensionless units, i.e., the pu system. To make the swing equation (47) dimensionless, we divide it by the rated power of the machine PR and define the quantities
The inertia constant H measures the ratio of stored kinetic energy and output power of the machine when it operates under normal conditions. Typical values are of the order of 4−6s. We then obtain the swing equation
where the power is now expressed in the pu system. The motion of the generator is coupled to the rest of the grid via the exchanged real power Pel, which we discuss later.
2. Principles of synchronous machines
As a next step toward a full machine model, we have to examine how the mechanical rotation induces voltages and currents. A careful introduction was given in Chap. 3.3 of Machowski, Bialek, and Bumby (2008); we provide here only the essential results. Consider first an unloaded round-rotor generator such as that depicted in Fig. 2. The rotor carries a field coil, where a dc current If generates a magnetic field. The stator consists of three coils called A, B, and C in which the rotating field induces an ac voltage; cf. Sec. III.A.1. The magnetic flux in the A coil with mutual inductance Mf, ΨrA(t)=MfIfcos(ωt), induces the voltage
commonly referred to as the “air-gap” electromagnetic force (EMF). The EMFs in coils B and C are the same up to a phase shift of ±2π/3. If the generator is connected to a load or the grid, a current
flows through the coil. It is generally phase shifted with respect to the voltage. Here I_=Ime−iφ/√2 is a complex (“phasor”) quantity. The voltage drop at the coil then reads
An explicit modeling of synchronous machines, including magnetic fluxes and dc currents, is cumbersome and not always necessary for analyzing large grids with tens to thousands of machines. Often an aggregated description in terms of voltage and current phasors is sufficient. Notable exceptions include the modeling of reactive power limits for voltage collapse, the detailed simulation of systemwide transients, or the design of power system stabilizers. In the aggregated description, one defines the two phasor variables E_f=ωMfIf/√2 and E_r describing the air-gap EMF. Equation (52) then reads E_r=E_f−iXAI_. We can further take into account imperfections, losses, and leakage of magnetic flux by an additional impedance Rl+iXl between the air-gap EMF and the terminal voltage of the generator V_g to obtain
Hence, the basic relation between the rotor quantities, the air-gap EMF at the stator coils, and the terminal generator voltage can be represented by a simple equivalent circuit, as depicted in Fig. 3. This description is effective and compatible with electric circuit theory.
In a salient pole machine the rotor is asymmetric; see Fig. 2. The asymmetry can be taken into account approximately in an effective model that differs only slightly from the case of a round-rotor machine. We decompose all phasors into two contributions, in and out of phase with the magnetic flux of the rotor. The two components are denoted by the subscripts d and q referring to the direct and quadrature axes, respectively. To describe the asymmetry, we replace the effective reactance XA with XAd and XAq for the d and q components. Hence, Eq. (53) is replaced by the two-component equations
with Xq=XAq+Xl and Xd=XAd+Xl. Figure 3 shows the equivalent diagram of the synchronous machine. Recombining the components yields the phasors of current and terminal voltage as follows:
Why have we introduced the EMFs Eq and Ep in addition to Ef? The quantities Ed,q characterize the physical source of the induced voltage (the magnetic flux of the field coil), while Ef quantifies the current in the field coil. The simple relations Eq=Ef and Ed=0 hold only during the steady operation of the machine. After a disturbance, the magnetic flux changes and so do the effective EMFs Eq,d. Hence, Eq,d become dynamic state variables for which we introduce the equations of motion in the next section. In contrast, Ef is a fixed system parameter or, if we include excitation control of the synchronous machine, a control variable.
3. Dynamics of synchronous machines
We now proceed to model the transient dynamics of a synchronous machine. A detailed derivation is available in textbooks on electric machines; see Chap. 11.1 of Machowski, Bialek, and Bumby (2008). Here we quote only the result and provide motivation for the form of the resulting equations.
Recall that a machine in the steady state was modeled by an EMF E_=E_q+E_d connected to the terminal ends via the reactances Xq and Xd, respectively. A similar approach can be used to describe the transient dynamics of a synchronous machine. The EMFs become time dependent, which is commonly indicated by a dash. Similarly, we must assign a different, transient value to the reactances. We stress that these are only effective quantities, and they yield a highly condensed description of the dynamics. A full dynamical model of a synchronous machine must be based on the magnetic fluxes in the different coils, the voltages, and the currents, and the model’s complexity increases quickly.
We first consider the EMF dynamics in the q axis, which quantifies the flux in the field coil (more precisely the field flux linkage). At steady state the flux is proportional to the voltage exciting the field coil, which leads to Eq=Ef. In the transient regime two additional terms have to be taken into account when one considers the voltage in the field coil. First, the magnitude of flux changes leads to an inductive term proportional to ˙E′q. Second, there is a term due to the coupling to the stator. A careful accounting of these terms yields the following equation (Machowski, Bialek, and Bumby, 2008):
The following similar relation can be found for the d-axis EMF:
We note that these two equations cover only the transient dynamics after a disturbance. Further models have been introduced to describe the short “subtransient” dynamics. We omit these models, as they typically yield only a small improvement (Stott, 1979; Weckesser, Jóhannsson, and Østergaard, 2013). Note, however, that higher-order models explicitly include damper winding, which we had to include phenomenologically in the swing equation.
4. Synopsis: A hierarchy of dynamical models
Combining the results from Sec. III.C.3 with the swing equation, we obtain a hierarchy of models; see Chap. 11.1 of Machowski, Bialek, and Bumby (2008). In the two-axis or fourth-order model the generator is described by four state variables, the transient EMFs E′q and E′d, the mechanical phase angle δ relative to the grid, and its derivative ω, which evolve according to
The one-axis or third-order model neglects the dynamics of E′d (Schmietendorf et al., 2014). Often it is set to zero ( E′d=0), as in the elementary steady-state model discussed in Sec. III.C.2. Hence, the model reduces to three state variables per generator and the last equation in Eq. (57) is omitted.
The second-order or classical model is widely used in the analysis of power system dynamics. It describes the mechanical motion of the generator at constant EMFs:
To study the resulting dynamics of the state variables, we have to specify the remaining quantities in the equations of motion. The electric quantities ˜Pel, Iq, and Id are time dependent. In fact, they are related to the state variables of all connected generators via algebraic equations describing the grid, as discussed later. The quantities Ef and ˜Pmech are set by the control system of the generator, more precisely, the exciter and the governor (Machowski, Bialek, and Bumby, 2008). If the action of the control system is ignored, they become fixed system parameters. The remaining quantities, in particular, the time constants T and the reactances X, are regarded as constant system parameters.
We note that detailed models can also include subtransient effects yielding a sixth-order model (Machowski, Bialek, and Bumby, 2008). Furthermore, dynamical models can be reformulated as adaptive networks (Berner, Yanchuk, and Schöll, 2021) or port Hamiltonian systems (Fiaz et al., 2013; Mehrmann et al., 2018), which enables further insights and generalizations.
5. A single generator coupled to an infinite busbar
First basic insights into power system stability can be obtained by focusing on the dynamics of a single generator; see Chaps. 5.3–5.5 of Machowski, Bialek, and Bumby (2008). To this end, we assume that the generator is connected to a large system and that the influence of the generator on that system is negligible. Such a setup is commonly referred to as an “infinite busbar” with a fixed voltage Vs. We carry out the calculation in the reference frame of the rotor and recall that the angle between the rotor and the stator frame is denoted by δ. Hence, the system voltages at the infinite busbar is written as
The generator is connected to the infinite busbar through a set of grid elements that are modeled as series impedances. There may first be internal losses or a leakage in the machine that is described by the impedance Rl+iXl. There is then typically a step-up transformer that we describe by the impedance RT+iXT, recalling that the pu system takes care of voltage levels automatically. Finally, there are transmission elements with effective impedance RS+iXS to the system. Typically, resistances can be neglected since Ohmic losses as well as shunt impedances are small. The machine EMFs and the system voltages are then related by
where we introduced the shorthand x′d=X′d+Xl+XT+Xs and x′q=X′q+Xl+XT+Xs. When one neglects Ohmic losses, the air-gap power acting on the machine equals the real power at the system terminal end
We thus found an explicit expression for the quantities ˜Pel, Iq, and Id in terms of the generator state variables, closing the equations of motion of Sec. III.C.4.
In the widely used classical model of power grid dynamics, one neglects transient saliency and thereby the asymmetry of the rotor geometry in a salient pole machine, assuming that x′d≈x′q=x′. This assumption simplifies the calculations, as we do not have to separate the q and d axes explicitly but can work with phasors. The relation of voltages and currents simplifies to
where the real parts correspond to the q components and imaginary parts to the d components. Writing E_′=E′e−iϕ and V_s=Vse−iδ the real power simplifies to
Finally, we recall that in the second-order model the dynamics of the rotor EMFs is fully neglected. Hence, E′ is set to a constant that is chosen to correspond to the steady-state value when the generator provides the power P+iQ. We then have to satisfy the relation
which can be solved for E2 with the result
6. Ohmic loads and the Kron reduction
To understand the operation and stability of extended power grids, including collective effects, we extend the model of Sec. III.C.5 to coupled generators. Recall that we have to specify the currents and the real power in terms of the internal EMFs in order to close the equations of motion summarized in Sec. III.C.4. To this end, we have to model the loads and we have to describe the grid explicitly. We formulate Kirchhoff’s equation in a fixed network frame of reference and work out the transformation to the reference frames of the rotating machines explicitly. For simplicity we neglect transient saliency here, setting x′d=x′q=x′. The derivation follows the presentation by Nishikawa and Motter (2015); a mathematical treatment was provided by Dörfler and Bullo (2013a).
Assume that the generator or active nodes are labeled as j∈{1,…,Ng} and the load or passive nodes are labeled as j∈{Ng+1,…,Ng+Nl}. For the generator nodes, we have already discussed the relation between internal EMFs E_′ and the voltage at the system terminal end V_ in Sec. III.C.5. The current injection from internal to terminal end at a node j is given by
Load nodes with a given power demand S_j=Pj+iQj are modeled by a fixed admittance to the ground ygroundj=S_∗j/V20. If the voltage magnitude Vj equals the reference voltage V0, the load nodes then consume the power V_jI_in*j=V0yground∗jV0=Pj+iQj, as desired.
For each node j∈{1,…,Ng+Nl} we then evaluate Kirchhoff’s current law: The current inflow must equal the current flowing to the other nodes in the grid:
We now collect all these linear relations in vectorial form, defining the vectors and matrices
Equation (66) and Kirchhoff’s current laws for the generator and load nodes are thus condensed into the form
where Ygg, Ygl, Ylg, and Yll are the respective partitions of the nodal admittance matrix (14). We can now gradually eliminate the network voltages, using the third row to eliminate V_l and then the second row to eliminate the V_g. This procedure is referred to as Kron reduction, a mathematical accounting of this procedure was discussed by Dörfler and Bullo (2013a). We finally obtain
where we defined the shorthand Y′=Ygg−Ygl(Yll+Ygrll)−1Ylg. Equation (70) has the same structure as Eq. (15), but the reduced susceptibility matrix Yeff represents an effective, not physical network. In summary, Eq. (70) directly yields the currents and the real power acting on the machines in terms of the internal EMFs as follows:
These quantities are given in the network frame of reference as previously stated. To evaluate them for the dynamical models summarized in Sec. III.C.4, we transfer back to the rotating frame of reference of each machine via
Solving for the currents in the machine frame yields
and the real power acting on the jth machine reads
This expression simplifies considerably for the particularly important case of a second-order model. We then have E′qj=Ej=const and E′dj≡0. Writing Yeffjk=|Yeffjk|ei(γjk+π/2) we obtain
7. The structure-preserving model
In Sec. III.C.6 all loads were modeled as constant admittances to the ground and subsequently eliminated. Bergen and Hill (1981) introduced an alternative model for power grid frequency dynamics, keeping the load nodes and hence the full network structure. Generator nodes are described as in the swing equations, setting the magnitude of the EMF to a constant. The model of load nodes is based on the observation that the real power consumption typically increases with frequency, which is written in the form
We thus obtain an equation of motion for load nodes that is equivalent to the swing equation with a vanishing inertia Hj. To close the model, one then has to express the real power injections Pelj at the generator and load nodes in terms of the voltage phase angles, recalling that the voltage magnitudes are assumed to be constant. This procedure, which includes the elimination of the generator terminal nodes, was described in detail by Nishikawa and Motter (2015). The great advantage of this approach is that the coupling in the resulting equations of motions reflects the true network structure.
8. Dynamics of power electronic inverters for renewable sources
Renewable power sources are generally not equivalent to synchronous machines. All photovoltaic power sources and virtually all wind turbines are connected to the grid via power electronic inverters. The transition to renewable energy supply thus fundamentally changes the dynamics of power grids. In particular, the decrease of inertia provided by conventional synchronous machines is a major challenge for grid stability (Milano et al., 2018).
One generally distinguishes between two modes of operation of power electronic inverters. A grid-following inverter provides a given amount of electric power adjusting to the voltage and frequency provided by the grid. Hence, the operation relies on other resources capable of providing a stable voltage and frequency. In contrast, grid-forming inverters regulate the voltage and frequency to specific set points similar to a synchronous machine. Mixed modes are also possible in principle, regulating either voltage or frequency and following the other.
The development and analysis of new types of grid-forming inverters and inverter-based grids is an active field of research. We thus introduce one important class, the droop-controlled inverter, following Schiffer et al. (2014). The basic state variables of such an inverter are the EMF magnitude Ej and the phase angle δj, where j labels the different generators in the grid. The control system adjusts these state values to maintain the predefined set values of power and frequency. To this end, the control system measures the real and reactive power exchanged with the grid and compares it to predefined set values. In a simple proportional, or droop control, scheme, the frequency control is proportional to the active power deviation and the voltage control is proportional to the reactive power deviation such that we obtain the equations of motion
where the superscript “mes” indicates the measured values of real and reactive power and the superscript d stands for desired values. Naturally, the desired frequency ωd is unique across the grid, whereas the desired voltages Edj may differ. The parameters κactj and κrecj are the droop gains for the frequency and voltage, respectively. The measurements are not instantaneous and are characterized by a low-pass filter such that Tj˙Pmesj=−Pmesj+Pelj and Tj˙Qmesj=−Qmesj+Qelj. The recovery time of voltage dynamics is typically much lower than the time constant Tj of the low-pass filter. We thus simplify the model by setting TVj=0 and obtain (Schiffer et al., 2014)
The real and reactive power exchanged with the grid ( Pelj and Qelj) depend on the state of all elements in the grid. To close the equations of motion, one thus has to specify the network equations. Modeling all loads by constant impedances to the ground as in Sec. III.C.6, we find that
where E_j=Ejeiδj; cf. Eq. (71).
9. Aggregated dynamical model
The classical model describes the dynamics of power systems as coupled second-order rotators, with each node representing one synchronous generator. Notably, similar equations of motion emerge on much coarser spatial scales (Zhang et al., 1997; You, Vittal, and Wang, 2004; Filatrella, Nielsen, and Pedersen, 2008; Chow, 2013), which emphasizes the generality of the coupled rotator model.
The dynamical model introduced by Filatrella, Nielsen, and Pedersen (2008) considered the aggregated dynamics of regions labeled by n∈{1,…,N}. If the internal coupling is sufficiently strong, we take the phase angle to be constant throughout the region. We again measure the phase and frequency relative to a frame rotating at the reference grid frequency ωR while writing ϕn(t)=ωRt+δn(t). The equations of motion for δn(t) are obtained from energy conservation. Each region stores kinetic energy in the rotation of all synchronous machines as follows:
where Jn is the aggregated moment of inertia. At each node a power Pinn is injected by generators or withdrawn by consumers and some energy will be dissipated at a rate Pdissn=ηn(˙ϕn)2. Assuming a constant voltage magnitude, the real power transmitted from a region n to a region m is determined by the phase difference
with Ohmic losses neglected. Knm is an aggregated quantity summing over all lines connecting the two regions n and m. The energy conservation for region n then gives
Noting that |˙δ|≪ωR in the vicinity of the normal operation, we simplify the expression for the change of the kinetic energy dEkinn/dt≈ωRJn¨δn and the dissipated power Pdissn≈ηω2R+2ηnωR˙δn. Substituting these results into the energy conservation law (83) then directly yields the equations of motion (Rohden, 2012)
where Peffn=Pinn−ω2Rηn. Equation (84) has the same structure as the classical second-order model but describes the dynamics on a coarser level. Equation (84) is based on elementary energetic arguments, but similar results can also be obtained via model reduction (Zhang et al., 1997; You, Vittal, and Wang, 2004; Chow, 2013). This emphasizes the generic nature of coupled rotator models in network science.
D. Load-frequency control
A power grid by itself does not store energy, so the generated power must match the demand and the losses at all times. A hierarchy of control systems exists to maintain this balance, which we introduce starting with the swing equation (47). As before, we use a rotating frame of reference such that ωi=˙δi denotes the deviation from the nominal grid frequency ωR,
In normal operation, all machines i∈{1,…,N} in a grid run in synchrony at the same frequency, a fact that we discuss further in Sec. VI. We thus focus on the bulk frequency and define
In many studies, the damping is assumed to be proportional to the inertia of a machines; cf. the Supplemental Material of Motter et al. (2013). Setting Di=ηJi we obtain
where ΔP=∑iPmechi−Peli is the power balance in the entire grid. We thus see that the power balance directly drives the dynamics of the bulk: A scarcity of generation leads to a decrease, while an overgeneration leads to an increase of the frequency. As the frequency can be measured easily, it is used to control the generation.
Different control mechanisms are distinguished according to the different timescales that they act upon and the different purposes that they serve; see Fig. 4.
(i) The momentary reserve is provided by the inertia of synchronous machines. The higher the inertia ¯J, the slower the grid frequency reacts to a power imbalance. This reaction rate is referred to as the rate of change of frequency (RoCoF) in the engineering literature. Notably renewable power sources with power electronic inverters have no intrinsic inertia, such that frequency stability is an important challenge in the energy transition (Milano et al., 2018).
(ii) Next, the power imbalance is reduced by increasing or decreasing the generation of specific control power plants or, more recently, battery electric storage systems (Fleer and Stenzel, 2016). According to the guidelines of the European Network of Transmission System Operators for Electricity (ENTSO-E) (UCTE Operations Handbook, 2009), primary control is activated as soon as the frequency leaves a small dead band around ωR and must be fully available within 30 s. Up to the dead band, the power is adapted proportionally to the frequency deviation as follows:
where the parameter κi describes the sensitivity of primary control for the grid area i. Primary control stabilizes the frequency but does not restore the nominal grid frequency, which can be seen by adding primary control to the equation of motion (87). For a constant power imbalance ΔP and a single area we find the fixed point
The damping terms ∝Di in Eq. (85) affect frequency dynamics similar to the primary control ∝κi, which is often directly included in the swing equation; cf. Sec. III.C.1.
(iii) Secondary control restores the nominal grid frequency ωR and reduces unscheduled power flows between different grid areas on a timescale of minutes. Control power plants adapt their generation according to the following proportional-integral (PI) law (Böttcher et al., 2020):
with κP and κI tunable gain factors. The local area-control error Gi is a measure of the power that is missing in area i. It is determined by the difference between the expected primary control power and the deviations ΔFi of the power flow to neighboring control areas, which is given as Gi=κiωi−ΔFi. Note that while the PI controller is linear, the local area-control error Gi depends nonlinearly on the system state. The finite time delay τ in the control [Eq. (90)] results from the time required for the measurement of Gi, communication, and the adaptation of the control power generation. The secondary control requires measuring local area-control errors Gi with a cycle time of 2–5 s (UCTE Operations Handbook, 2009). The delay time τ can therefore be expected to be of similar magnitude but may vary in different control areas and with time (Böttcher et al., 2020).
(iv) Tertiary control power is invoked on even longer timescales. In Europe, the tertiary control power must be activated within at most 15 min. Tertiary control is partly manually operated and used mostly to restore the automatic control reserve.
Load-frequency control scheme according to ENTSO-E (UCTE Operations Handbook, 2009).
IV. Power Grid Topologies and Datasets
We previously introduced the various components of power grids, in particular, power lines and assorted producers and consumers at the nodes. To derive a model of a complete power grid, the missing ingredient is the topology, that is, how these components are connected. Besides studying concrete examples of real-world power grids, it is also desirable to obtain plausible synthetic networks to allow a more systematic study of the impact of topological features on dynamical properties.
A. General aspects
Large power grids, especially at the continental scale, operate at many different voltage levels. Long-range transmission is achieved through high voltage connections in order to minimize losses. In Central Europe, for example, the transmission grid operates at 220 and 380 kV. The regional distribution is then achieved by grids at successively lower voltages that are connected to the high voltage transmission grid at substations (transformers). Typical voltage levels in Europe range from high voltages (such as 60 and 110 kV) over medium voltages (3–30 kV) to low voltages (such as 230 and 380 V). Typical line parameters for transmission grids are summarized in Table I.
Typical parameters of overhead transmission lines at different voltage levels for the United Kingdom and the USA according to Kundur (1994) and Machowski, Bialek, and Bumby (2008) and for Germany according to Oeding and Oswald (2016). Ohmic resistance, susceptance, and charging capacity are proportional to the length of the line and are given per km.
Grids at these different voltage levels are constructed in different ways and with different trade-offs and goals in mind and thus exhibit different network topologies. A low voltage distribution grid is typically organized as a tree graph, which in this context is called a radial design. While loops might physically exist, there often is a switch interrupting them. We remark that only if a line fails is the loop closed to resume supply downstream of the failure.
High voltage transmission grids require uninterrupted operation even in cases where a line fails (the N−1 criterion). This implies that the network is meshed; i.e., it contains loops. An example of such a topology is shown in Fig. 5. Medium voltage grids often fall in an intermediate regime, with varying amounts of meshing.
Topology of a high voltage power grid. The visualization has been produced by Paul Cuffe using the methods described by Cuffe and Keane (2017) based on data of the NESTA archive. From Coffrin, Gordon, and Scott, 2014.
From a network science perspective a notable aspect of power grids is that they are geographically embedded. That is, nodes and lines have geographical locations and lengths. Line crossings are possible, but rare. These aspects are reflected in the connectivity or topology of the network, putting them in a different class than the ensembles typically studied in network science.
Further, as infrastructure has costs, power grids typically also contain as few lines as possible. As a consequence, power grids are typically sparse. The observed average degree for transmission grids falls between 2 and 5, and the degree distribution peaks between 2 and 3 and has an exponential tail (Solé et al., 2008; Wang, Scaglione, and Thomas, 2010).
B. Network ensembles and synthetic grid models
The statistical physics of networks often considers ensembles of networks to map out typical topological properties and to provide a reference class for actual real-world datasets. Power grids do not fully resemble any of these common network ensembles. In particular, they are not well captured by small-world networks (Cotilla-Sanchez et al., 2012). Although the clustering coefficient is similarly high, there is a much more pronounced connectivity at large scales. A Watts-Strogatz small-world network of such a low degree would likely end up being disconnected for moderately sized networks, and the nonlocal links in power grids tend to be much more clustered (Wang, Scaglione, and Thomas, 2010).
To capture the full topological structure of power grids, we need to define novel network ensembles: probability distributions on the space of networks that capture the notion of what looks like a power grid and what does not. This is typically done by specifying a random process that produces networks, either by rewiring links or by iteratively growing networks. Given such an embedded topology, typical power line parameters can then be used to derive the power grid’s properties (Table I).
Wang, Scaglione, and Thomas (2010) proposed a rewiring-based process with the aim of staying close in spirit to the construction of small-world networks while taking some specific identified properties of power grids into account. In contrast, growth-based models, starting with that of Schultz, Heitzig, and Kurths (2014b), grow the grid with random node placements and new connections. The model of Schultz, Heitzig, and Kurths (2014b) mimics the trade-off between global resilience and economy in the growth process and is able to recreate the exponential degree distribution. Stating these trade-offs explicitly allows one to study the trade-offs of dynamic and structural stability given by Plietzsch et al. (2016). A growth model for hierarchical networks was proposed and analyzed by Ódor and Hartmann (2018).
Soltan and Zussman (2016) considered node placement in addition to line generation. They considered average path length, clustering coefficient, the slope of node degree distributions, and the line length distributions and showed that their algorithm matched those in various North American power grids.
In these models the growth process makes use of topological and embedding information but does not consider the resulting energy flows in the growth stage. Birchfield et al. (2017) and Birchfield, Xu, and Overbye (2018) added consideration of power flows in the dc approximation and voltage profiles to the iterative growth of synthetic grids. Soltan, Loh, and Zussman (2019) fitted the growth process using a Gaussian mixture model. Espejo, Lumbreras, and Ramos (2019) presented a model that focuses instead on the historical plausibility of the growth process and economy versus robustness trade-offs.
Besides generating topologies from scratch, it is also possible to take an existing topology and vary it while keeping many aspects of it fixed. Such surrogate ensembles for spatial infrastructure networks were described by Wiedermann et al. (2016). A difficult open question is to find a fully satisfactory model for the cogeneration of a variety of networks (Hackl and Adey, 2019), or hierarchies of networks (Schultz et al., 2016).
C. Network datasets
The availability of real grid data is limited. Thus, a variety of test cases have been compiled for scientific use. Some are based on the available data of real-world grids. Some are synthetic ones trying to mimic the structural properties of real-world grids. In the following, we list resources that provide grid data without a claim of completeness.
• The classical IEEE test cases are probably the most heavily studied power grids of all, with some having been implemented decades ago. The datasets are partly synthetic and partly derived from real grids and can be obtained at the repository (Christie, 1999).
• Several online repositories have been created in recent years to provide larger, more diverse, and more recent test cases focusing on transmission grids (Birchfield, 2016; Pacific Northwest National Laboratory and National Rural Electric Cooperative Association, 2017; Farid et al., 2019) or distribution (“feeder”) grids (Schneider et al., 1991; Kavasseri and Ababei, 2021). In addition, several papers reviewed such repositories and the methodologies; see Coffrin, Gordon, and Scott (2014), Birchfield et al. (2017), and Schneider et al. (2018). These repositories again include both synthetic grids and approximate models of real grids, with a focus on North America. Large test grids with properties typical for European grids were provided by the Pegase project (Villella et al., 2012).
• The IEEE and Pegase test grids are included with others in the latest version of the software package matpower (Zimmerman, Murillo-Sanchez, and Thomas, 2011).
• Many test cases are restricted to static properties and do not include the parametrization necessary for dynamic simulations. To analyze the dynamics of these systems, some heuristic assumptions about present or future parameters are necessary. Realistic dynamic test cases mostly come in the form of small synthetic test cases developed by IEEE that can be found in repositories (Pacific Northwest National Laboratory and National Rural Electric Cooperative Association, 2017; Farid et al., 2019). ENTSO-E provides semisynthetic grid models mimicking the European grid (Semerow et al., 2015), but access is limited. A notable exception is the synthetic model of the northeastern U.S. published by Birchfield (2016), which includes a 25 000-node dynamic system, validated in the sense of the system of Xu, Birchfield, and Overbye (2018). A heuristic dynamical parametrization of the ENTSO-E–based system (Wiegmans, 2016) was recently given by Pagnier and Jacquod (2019c).
• Recently several initiatives started to map out power grids from publicly available data. A model of the German power grid was extracted from OpenStreetMaps (Medjroubi, Matke, and Kleinhans, 2015), and a model of the European grid was extracted from the ENTSO-E interactive grid map (Wiegmans, 2016). A variety of related datasets can be found on community websites (Open Energy Modelling Initiative, 2017; Open Power Systems Data, 2017). A collection of grid frequency time series was presented by Gorjão, Jumar et al. (2020).
V. Dynamics of Elementary Networks and Building Blocks
We introduce basic aspects of power system dynamics and stability for an elementary grid containing one transmission line. We first consider its static operation and its limitations before we turn to the dynamic stability.
A. Static solutions and voltage stability
We consider the steady state of an elementary circuit with one load node and one generator connected by a single transmission line to understand which factors limit the power transmission in different types of grids. The generator s is treated as a slack node, so Vs=1 and δs=0 are fixed; cf. Sec. III.A.3. Most load nodes draw the reactive power Qn at a fixed ratio to the real power, which is commonly specified in terms of the power factor cos(ϑ) defined via Pn=−√P2n+Q2ncos(ϑ) and Qn=−√P2n+Q2nsin(ϑ). We analyze the operation here as a function of the real power demand Pn<0.
To obtain the voltage magnitude Vn and phase angle δn at the load node, we then solve the ac load flow equations (25), which are rewritten using the addition theorems of the sine and cosine
where the angle γns is defined via
We can then separate the equation for the voltage magnitude using sin2+cos2=1, which yields
which is easily solved for V2n. The phase angle can then be obtained by solving one of Eqs. (91).
To understand the behavior of the system and the limitations to power transmission in more detail, we plot the voltage Vn as a function of the real power demand |Pn| in Fig. 6(a) for different values of the power factor. This plot is often referred to as the “nose curve” in power engineering (Machowski, Bialek, and Bumby, 2008). First, one observes that the quadratic equation (92) can have two solutions, where only the branch with the higher magnitude is relevant, as the other branch is unstable. Second, power transmission is generally limited. Physical solutions of Eq. (92) exist only if the real-power demand does not exceed the limit,
which depends on the reactive power. The limit as well as the voltage magnitude Vn is lower when the load draws reactive power [ Qn<0, tan(ϑ)>0]. However, if the load node supplies reactive power instead [ Qn>0, tan(ϑ)<0], the voltage at the load node can be even larger than Vs and the transmittable real power is strongly enhanced. Grid operators hence aim to keep tan(ϑ) small or even negative by supplying reactive power near the loads. This procedure is referred to as “reactive power compensation.”
Nose curve. (a) Voltage at the load node as a function of transmitted power for three different power factors: tan(ϑ)=+0.3,0,−0.3 (from left to right). (b) Voltage at the load node for tan(ϑ)=0 in the initial state (right) and after half of the connectivity is lost, i.e., Bns is halved (left). Parameters are Sbase=100 MW, Vs=1 pu, and Bns=10 pu. Solid (dashed) lines indicate stable (unstable) solutions of the load flow equations.
We further consider how a “voltage instability” can arise in elementary power grids. We assume that two parallel transmission lines exist, and that one of them is lost in a contingency situation. As a consequence, the effective susceptance Bns is halved and the nose curve contracts, as shown in Fig. 6(b). Further effects depend on the grid loading. If the real power demand |Pn| is not too high (such as the 200 MW value in the figure), a solution of the power flow equations still exists, albeit at a lower voltage. The consumers will thus experience a rapid drop in the voltage level. However, if the real power demand is too large (such as the 300 MW value in the figure), no static equilibrium solution exists anymore and a “voltage collapse” emerges in which voltages decline as governed by the power system dynamics. A comprehensive analysis of voltage stability in power grids was given by Van Cutsem and Vournas (2007).
B. Flow and voltage limits
Power transmission in current highest-voltage grids is typically limited by factors other than the ones discussed earlier. Most grids are constantly monitored and thoroughly regulated. Security limits have been formulated for current and voltage [see European Network of Transmission System Operators for Electricity (2004)], and emergency shutdowns are carried out when these limits are violated. Hence, the physical limits of power transmission are rarely met during normal operation.
In particular, extreme currents can lead to an Ohmic heating and eventually to a dangerous bending of overhead transmission lines due to thermal elongation. During the 2003 power outages in North America (U.S.-Canada Power System Outage Task Force, 2014) and Italy (Union for the Co-ordination of Transmission of Electricity, 2004) transmission lines hit trees, leading to a short-circuit fault. Hence, grid operators typically impose “thermal limits” for the currents and overloaded lines switch off automatically. Similarly, upper and lower limits for the voltage magnitude are imposed to guarantee power quality and avoid the danger of a complete voltage sack.
To see how these regulations limit power transmission in highest-voltage grids, we solve the load flow equations (91) as a function of Pn, assuming a fixed power factor set to a typical value of cos(ϑ)=0.95, with ϑ>0. We use typical parameters for a 380 kV transmission line in Western Europe with r/l=0.03 Ω/km and x/l=0.246 Ω/km and a thermal limiting current of Ith=2.58 kA (Oeding and Oswald, 2016), with l the length of the line. We use the notation of Sec. V.A but give all quantities in physical units for better accessibility.
Figure 7 shows the results for two characteristic cases. We first consider a transmission line of length l=100 km connecting the load node and a generator node, which has a fixed voltage magnitude Vg=400 kV, i.e., slightly above the nominal voltage level of 380 kV. The current
increases approximately linearly with the transmitted power. For |Pn|≥1.61 GW, the current exceeds the thermal line limit, which would lead to emergency measures, whereas the voltage magnitude stays well above the lower limit of Vlim=0.9×380 kV. This illustrates that in densely populated areas with short lines, it is mainly the thermal limit for the current that limits the operation of a transmission line. In fact, we find a good estimate for the maximum transmittable real power using
Voltage stability can be an issue if generators are located far from loads. Figure 7 shows the current and voltage at the load node for a longer line with l=200 km. The voltage limit Vlim is already hit for |Pn|≥1.01 GW, where the current is still well below the thermal limit.
Factors limiting real power transmission in ac power grids for two different characteristic cases (blue, short line with l=100 km; dashed red, long line with l=200 km). We consider a single transmission line connecting one power plant to a substation. Shown are the current |I|, the voltage magnitude |Vn|, and the voltage angle δn at the load node as a function of the real-power load |Pn|. The dotted line shows the limits Ith and Vlim for the current and the voltage.
Voltage limit violations can be handled to some extent by using reactive power compensation, as discussed in Sec. V.A. Decreasing the power angle ϑ at the load node allows for a higher power transmission. In contrast, the current limit represents a more severe limitation. A simple way to extend power transmission is dynamic line rating, where the thermal limiting current Ith is adapted to the ambient conditions (Douglass and Edris, 1996). On cold windy days, a higher current is acceptable, without the risk of overheating.
C. Dynamic stability
1. Local stability of a single generator
We now turn to the question of dynamical stability, starting with the most elementary setup: a single generator coupled to an infinite grid. The grid is assumed to be so large that it is not affected by generator dynamics, such that its voltage magnitude V and phase angle ϕ are fixed. We consider here the question of angular stability, and thus assume that the generator voltage is constant at the steady-state value E∘. If the phenomenon of voltage collapse occurs, voltages drop dynamically such that a more refined treatment is necessary.
Assuming a lossless connection of the machine to the grid with the susceptance x, the real power flow from the machine to the grid is given by Pel=Ksin(δ−ϕ) and the swing equation (49) reads
where K=VE∘/x is an effective coupling strength. In normal operation the machine rotates with the frequency ωR and a fixed phase difference to the grid to supply a constant real power. Hence, we are interested in the fixed points of Eq. (96) given by
One finds that the + solution is linearly stable, while the other branch is unstable. The solutions vanish in a saddle node on a circle bifurcation when Pmech=Pcrit=K. Hence, the power that the machine can receive and transmit to the grid is limited to a maximum value of Pcrit.
The nature of this bifurcation becomes more obvious with a mechanical analog. Consider a point particle with an effective mass meff=2H/ωR and a friction coefficient η=D moving in a tilted washboard potential
where δ denotes the particle position. Newton’s equation for the particle meff¨δ=−η˙δ−∂Veff/∂δ is then equivalent to the swing equation (96). The tilted washboard potential is shown in Fig. 8 for different values of tilting Pmech. For low tilting, stable (unstable) fixed points exist in the local minima (maxima) of the potential (green and red dots). As the tilting increases, the minima and maxima approach each other and finally vanish for Pmech>Pcrit. The tilted washboard potential has been studied in great detail in statistical physics, particularly in the presence of noise; see Risken (1996) and references therein.
The tilted washboard potential [Eq. (98)] as a mechanical analog to the dynamics of a single generator in the second-order model for (a) Pmech=0.1 K, (b) Pmech=0.5 K, (c) Pmech=0.9 K, (d) Pmech=1.1 K. The potential Veff is given in units of K. Green disks, stable fixed points; red crosses, saddle points.
2. Bifurcations
To obtain a comprehensive understanding of the dynamics of a single generator, we also analyze its global stability properties. Given an initial state δ(0),˙δ(0), what happens for long times? Does the generator relax to the previously discussed stable fixed point?
A global phase-space portrait of the single machine system is shown in Fig. 9 for different values of Pmech/K. One observes that in most cases a limit cycle exists that corresponds to a desynchronized generator. The generator cannot exchange real power with the grid in this state, as Pel=Ksin(δ−ϕ) is oscillating and averages out to zero. The mechanical input power accelerates the generator until this input is balanced by the damping such that it has a higher frequency than the grid ( ˙δ>0).
Basin of attraction (red) of the stable fixed point of Eq. (96) in the δ- ˙δ plane. Parameters are 2H/ωR=1, D=0.3, K=8, and (a) Pmech=1, (b) Pmech=2, (c) Pmech=4, and (d) Pmech=9. The white dots are the stable and unstable fix points. The black line gives the limit cycle.
We thus find three parameter regions with different global stability properties; see Figs. 9 and 10. For Pmech/K>1 no fixed point exists and the system will always converge to the limit cycle. If Pmech/K is small or the damping η is large, the system is globally stable: It converges to the attractive fixed point for almost all initial states. In the remaining part of the parameter space, the fixed points and the limit cycle coexist such that the long-term dynamics depends on the initial state.
Stability map of a single generator coupled to an infinite grid in the second-order model. Three different parameter regions exist: a globally stable fixed point, a globally attractive limit cycle, and a region of coexistence of a fixed point and a limit cycle. Bifurcations between these regions are shown as colored lines. The red dashed line is the approximation (102) for low friction. See Risken (1996) and Manik et al. (2014) for a more detailed analysis.
An analytical approximation for the border between the globally stable and the coexistence regime can be obtained in the low-friction limit using Lyapunov’s second method (Parks, 1992; Risken, 1996; Manik et al., 2014). We define an energy functional E=meff(˙δ)2/2−Kcos(δ), setting ϕ=0 for simplicity. Using the equation of motion for δ, we obtain
If E decreases on average for all initial conditions, the system is in the globally stable regime. The condition for the border between the globally stable and the coexistence regime is therefore obtained by setting ¯¯¯¯¯¯¯¯¯¯dE/dtT=0, where the bar denotes the average over one period T. Evaluating this condition yields
We can now calculate ¯˙δT=2π/T and obtain
At the bifurcation, when a globally stable fixed point loses stability, there is a trajectory that satisfies ˙δ=0 at each successive peak of the potential landscape. This trajectory has Epeak=K, which can be assumed to be constant in the low-friction limit. Replacing E(δ,˙δ)≈K, the integral in Eq. (101) can be evaluated in closed form and the border between the globally stable and the coexistence region [Eq. (100)] is given by
Equation (102) matches the results from a direct numerical evaluation for small values of η, as shown in Fig. 10. For large values of η it slightly overestimates the globally stable parameter region.
3. Probabilistic stability measures
Elementary systems allow for a comprehensive understanding of the bifurcation structure and global stability. However, this approach is no longer feasible for high-dimensional networked systems, such that more versatile methods are needed. Recently, probabilistic approaches have attracted strong interest (Menck et al., 2013). For instance, one can quantify the stability of a dynamical system by the probability of returning to the desired fixed point after a random perturbation. The set of all initial states guaranteeing convergence, the basin of attraction, is shown in Fig. 9 for a single generator and different system parameters. While the geometry of the basin can be rather intricate in higher dimensions, its size is easy to measure by randomly drawing a sufficient number of initial states from a predefined set A and numerically simulating their dynamics. It is then the same as the following Monte Carlo integration of the indicator function 1B on the basin of attraction:
Different aspects of stability can be understood by choosing an appropriate set A. For instance, restricting perturbations to a single node reveals weak spots in the grid (Menck et al., 2014).
The process of sampling such a probability is a Bernoulli process. Its sampling uncertainty decreases with the number of samples N, which is asymptotically given by the binomial confidence interval
although better estimates are possible for small values of N (Agresti and Coull, 1998). Notably the number of samples needed does not increase with the dimension of the problem, whereas the individual samples are typically more expensive to obtain.
The asymptotic convergence to a fixed point is an important aspect of stability, but it is not sufficient in the context of power grids. One can generalize the current approach, requiring that a trajectory does not exceed certain operational limits or, equivalently, that it does not not leave a “survival” region S. This leads to the following definition of the survivability (Hellmann et al., 2016):
Survivability requires the specification of a survival region S, as shown in Fig. 11 for a single generator, assuming that the survival region is given by |˙δ|≤±5 rad/s. This assumption is useful in the study of power grids, as generators have many protection circuits that would switch them off to avoid mechanical damage. Probabilistic measures have also been generalized to cover the case of repeated perturbations and stochastic systems (Schultz et al., 2018; Lindner and Hellmann, 2019). Furthermore, the uniform sampling from the set A can be replaced by more general distributions of initial states.
Survivability region (red) of the stable fixed point of Eq. (96), with the survival region given by the area between the lines ˙δ=±5 rad/s. The basin of attraction of Fig. 9 is shaded in orange. The parameters are 2H/ωR=1, D=0.3, K=8, and (a) P=1, (b) P=2, (c) P=4, and (d) P=9. The white dots are the stable and unstable fixed points. The black line gives the limit cycle.
4. Stochastic stability
Thus far we have analyzed the stability of the swing equation for fixed system parameters. But real-world systems will typically face disturbances and noise. How does this affect the stability of the system?
We now analyze how the stability of the swing equation (96) is affected if the mechanical input power of a machine fluctuates in time as
Notably this problem can be solved almost fully analytically in the case of white noise ξ(t) using Kramers’s escape rate theory (Risken, 1996). Using the analogy to a particle moving in a tilted washboard potential [Eq. (98)], the main question is as follows: Can the particle overcome the potential barrier and thus escape the vicinity of the stable fixed point as illustrated in Figs. 12(a) and 12(b)? What is the probability of such a desynchronization event?
Desynchronization of the swing equation due to a fluctuating power input. (a),(b) When the input power Pmech fluctuates, the generator can lose synchrony to the grid after an escape time τ. (c) Kramers’s escape rate theory predicts the escape process for Gaussian white noise. The theoretical prediction [Eq. (107) (black lines)] accurately predicts the mean escape times ¯τ obtained from direct numerical simulations averaged over 2000 sample processes (symbols). Results are shown for an extremely loaded connection with K=1 and P=0.95 and a damping coefficient of D=8×10−5 s−1×ωR (intermediate damping). From Schäfer et al., 2017.
Kramers’s theory gives the average time ¯τ until a particle escapes (van Kampen, 2007; Schäfer et al., 2017) as follows:
where σ2 is the variance of the Gaussian white noise, η=D is the effective friction, and 2ϱ=η+√η2+(8H/ωR)|V′′eff(δmax)| for intermediate damping. The quantities V′′(δmin,max) are the second derivatives of the potential evaluated at the local minimum and maximum δmin,max, respectively, and
is the potential barrier that the particle has to overcome. Numerical simulations show an excellent agreement with this prediction, as indicated in Fig. 12(c).
Kramers’s formula (107) reveals how the system parameters affect the stability. The main factors are noise strength σ, damping parameter D, and potential barrier ΔVeff, which enter Kramers’s formula exponentially. In contrast, ¯τ depends less sensitively on the inertia H, which enters only algebraically.
The barrier is determined by the effective tilting of the potential, i.e., the loading of the grid, and vanishes as ¯P→Pcrit. Figure 12 demonstrates the escape process for an extreme load, whereas ¯τ is longer by orders of magnitude for realistic loading levels. However, ¯τ can be drastically reduced if the noise is no longer Gaussian but instead intermittent, which is typical for wind power systems (Schmietendorf, Peinke, and Kamps, 2017). In larger power grids, different escape routes can exist, corresponding to different parts of the grid loosing synchrony (Schäfer et al., 2017). In Sec. VIII.D we reconsider the stability of the swing equation under several generalizations.
VI. Synchronization and Steady States in Complex Grids
A. The need for synchrony
The stable steady operation of a power grid requires perfect phase locking of all synchronous machines throughout the grid (Rohden et al., 2012; Dörfler, Chertkov, and Bullo, 2013; Motter et al., 2013). That is, all machines in a connected network have to run at the same frequency with a well-defined phase difference between them; otherwise, no steady power flow is possible. From Eq. (73) we find that if the nodes operate at their nominal voltage amplitude Ej=1, the real power flow between two nodes i and j is determined by their phase difference
where Yeffij=Kijei(γij+π/2) is the effective nodal admittance and the term sin(γij) comes from the diagonal of Yeff. The same expression is used for the aggregated model introduced in Sec. III.C.9 with γij=0.
Under fault conditions, a fast transient instability can occur when a group of machines accelerates relative to the others. The relative phases δi,j(t) between the groups are then no longer locked and there is no steady power flow Pi→j. If the fault clears, another transient may either diverge or converge. Transient instabilities on timescales of cycles typically lead to the shutdown (“tripping”) of both transmission lines and generators.
Furthermore, violations of perfect phase locking are observed during “interarea oscillations,” where electric power oscillates across the grid at low frequencies of 0.1–10 Hz (Klein, Rogers, and Kundur, 1991; Rogers, 2012). Inter-area oscillations repeatedly occur after disturbances such as the loss of a generator and are typically damped out in minutes. However, oscillations may also grow in exceptional cases, eventually leading to a system split or blackout. For instance, a cascade of failures led to the tripping of several transmission lines and the loss of the stable steady state during the 2006 European power outage (Union for the Coordination of Transmission of Electricity, 2007). Oscillations grew rapidly and the grid finally separated into three mutually asynchronous areas.
We are thus led to a fundamental question of power system stability: Do the equations of motion admit a stable solution where all nodes are perfectly phase locked,
given the network topology and the power injections? Synchrony is about phase dynamics; therefore, it is natural to focus on the phases as dynamic variables in the following. To simplify the analysis, we neglect voltage dynamics for the time being and assume that the phase dynamics at the nodes is well described by the swing equation (49). As previously discussed, this equation is generic and arises similarly for single synchronous machines, grid-forming inverters with power control, and aggregated dynamical models. Absorbing the various factors in front of ¨δi and ˙δi into reduced inertia Ji and damping Di coefficients, the central equation in question then reads
The general problem of synchronization and dynamical stability with voltage dynamics is more involved, and we present an outlook at the end of this section.
Versions of this model have also been studied in the theoretical physics literature under the name of the second-order Kuramoto model; cf. Chap. 7 of Rodrigues et al. (2016). However, note that different aspects of synchronization have traditionally attracted more attention in theoretical physics. Starting from the seminal work of Kuramoto (1975) many researchers have investigated the transition from incoherence to partial synchronization, where some nodes lock their frequencies. This is not sufficient for power systems. Furthermore, many essential results on partial synchronization have been obtained in a large- N mean-field limit. A generalized large- N limit for complex networks and the problem of full synchronization were treated only recently by Kuehn and Throm (2019).
B. Linear stability of the phase dynamics
We begin by assuming that a solution to the static equations is given. We then study whether its phase dynamics is stable to small perturbations. Consider a small perturbation around the synchronous fixed point δi(t)=Ωt+δ∘i+αi(t). To linear order, perturbations evolve according to
with the matrix Λ∈RN×N defined as
The fixed point is linearly stable if all perturbations α are damped. We note that power grid synchrony is not affected if all angles δi are shifted by a constant. Thus, this mode is excluded from the stability analysis. One can show that stability is primarily determined by Λ: The synchronous state is linearly stable if all eigenvalues of Λ are positive except for the one corresponding to a global shift of all phases. This criterion becomes especially simple if line losses can be neglected. γij=0, and Λ is a weighted, symmetric Laplacian matrix; cf. Sec. II. Stability is guaranteed as long as the “weights” Kijcos(δ∘i−δ∘j) are positive for all links (i,j); i.e., if no line is overloaded,
This criterion is sufficient, but generally not necessary. It can be viewed as a generalization of the stability criterion for a single generator discussed in Sec. V.C.1.
The situation becomes more involved when line losses cannot be neglected. The condition |δ∘i−δ∘j−γij|<π/2 is then not sufficient and more subtle necessary conditions are derived (Skar, 1980). Further effective criteria for the definiteness of signed Laplacians were given by Song, Hill, and Liu (2015) and Chen et al. (2016). Lyapunov exponents and vectors for this system were studied by Bosetti and Khan (2018). Note that phase differences in transmission grids are typically much smaller. Other factors limit power transmission, at least for short transmission lines; see Fig. 7.
Stability can be lost if a parameter of the network is varied, e.g., if power injections and grid loads increase or transmission lines fail (Coletta and Jacquod, 2020). Eventually, the second Laplacian eigenvalue tends to zero and the fixed point is lost in an inverse saddle-node bifurcation on a circle (Manik et al., 2014). Notably bifurcation sets and linear stability properties become much more intriguing if voltage dynamics is included (Schmietendorf et al., 2014; Ma et al., 2016; Sharafutdinov et al., 2018) or in inverter-based grids (Schiffer et al., 2014). An overview of different aspects of dynamical stability was provided by Gajduk, Todorovski, and Kocarev (2014).
C. Existence, multiplicity, and properties of synchronous states
We now return to the question of whether a stable synchronous state exists for a given power grid. We first review some important properties of synchronous states and then review a theoretical approach to systematically compute and classify all such states.
1. Properties of synchronous states
The proper operation of a power grid requires perfect phase locking according to Eq. (110). Substituting into the equations of motions (111) yields the following fundamental condition for the phases and active power flows in the synchronous state of a power grid:
We note that the flows are asymmetric for γ≠0 due to Ohmic losses. In fact, one can express the losses at a transmission line (i,j) as
Losses vanish when γij=0 and when there is no power flowing on the line, i.e., when δ∘i=δ∘j.
Synchrony does not necessarily imply that the grid operates exactly at the reference frequency of 50 or 60 Hz. In fact, we can determine the equilibrium frequency directly from Eq. (115). Summing over all nodes i∈{1,…,N} and solving for Ω yield
We recall that the effective parameters Di include both damping and the action of primary load-frequency control, such that the current result is equivalent to Eq. (89). The deviation from the reference frequency Ω is given by the ratio of the power imbalance, including losses, and the cumulative primary control strength. We return to the analysis of the bulk frequency in Sec. VIII.A.
2. Phase cohesiveness and necessary conditions
We are particularly interested in synchronous states where angle differences are small; cf. the discussion in Sec. V.A. Small angle differences are also essential for many analytic results on the stability of the swing equation; cf. the linear stability condition (114). To state this aspect more precisely we define the cohesiveness as follows: A state is called phase cohesive if all phases δ∘i lie in an arc of length ζ∈[0,π/2) denoted as A(ζ).
The definition can be applied to deriving necessary conditions for synchronization. For simplicity, assume that the power in the grid is balanced and that we have no losses ∑iPi=0, γij=0. Hence, we have Ω=0 in Eq. (110) and the condition for phase locking in the second-order power grid model (111) reads
For a phase cohesive sync state, the magnitude of the sum can be bounded from above by sin(ζ)∑nKin such that Eq. (117) can be satisfied only if, for all i∈{1,…,N},
Furthermore, we can consider any two nodes i and j of the grid and evaluate the difference of the conditions (117). We can again bound the sums on the right-hand side which yields the following condition for all i,j∈{1,…,N}:
The two conditions are necessary for a cohesive sync state to exist, but they are not sufficient by any means. Nevertheless, they reveal two important facts: (i) Synchronization requires a sufficient connectivity of the grid, while a strong divergence of the Pi impedes synchronization. (ii) If the phase divergence ζ is reduced, we have to increase the coupling or decrease the divergence of Pi. Further generalizations and refinements were discussed by Chopra and Spong (2009), Ainsworth and Grijalva (2013), and Dörfler, Chertkov, and Bullo (2013).
3. Existence of solutions and multistability
We now turn to the question of whether solutions to the power flow equations exist. Notably power grids can admit multiple stable and unstable states of operation, even when power injections are fixed. Sudden changes of the power injections or the grid topology can trigger transitions between these states, with drastic consequences for grid operation and stability. Circulating power flows can emerge, which are generally undesired, as they increase line loads and Ohmic losses. A classic example of such circulating power flows is the Lake Erie loop in North America (Coletta et al., 2016; Jafarpour et al., 2019). If the system finds itself in a situation in which no stable synchronous state exists, system collapse is an inevitable consequence. We begin with the simplest case of a power grid without Ohmic losses such that the fixed points are determined by Eq. (117). The central questions now are as follows: When does this equation have multiple solutions that correspond to a dynamically stable power flow, and why?
We start with an elementary example: a simple ring network where Pi=0 for all nodes i∈{1,…,N} and homogeneous line parameters K (Schröder, Timme, and Witthaut, 2017). One can easily verify that the fixed points of this network are given by δ∗i=2πim/N, with a parameter m∈{−N/2,−N/2+1,…,+N/2}. Fixed points with m∈(−N/4,+N/4) satisfy Eq. (114) and are guaranteed to be dynamically stable; all other fixed points are potentially unstable. This simple example already suggests the fundamental features of multiplicity of sync states. Power flow solutions in lossless grids can differ in cycle flows. In fact, each fixed point corresponds to a flow f=Ksin(2πm/N) around the cycle. The number of stable fixed points increases with the size of the loop and is limited to one if N≤4.
Using these insights, we introduce a general method to construct all stable fixed points of a generic lossless power grid following Manik, Timme, and Witthaut (2017). The main idea is to shift the focus from nodes to the edges and cycles of the network. Hence, we define a vector of flows on the network’s edges as
where E denotes the node-edge incidence matrix (4) and the sine function is taken elementwise. The steady-state conditions (117) then become
Fixed points are now obtained via a stepwise procedure: We first construct all solutions of Eq. (121) and then reject those that are incompatible with Eq. (120).
The first step is based on the following observation: The kernel of the incidence matrix E corresponds exactly to cycle flows in the network. Hence, the general solution of Eq. (121) can be written as
where F(s)∈RL is a specific solution, f∈RL−N+1 gives the strength of the cycle flows, and the matrix C is the cycle-edge incidence matrix (5). For each flow vector F, we can construct the associated nodal phases as follows. Start with the slack node s and set δ∘s=0. Proceed to a neighboring node j. If one assumes that the connecting edge e^=(j,s) is oriented from node s to node j, the phase value reads δ∘j=δ∘s+Δe, where the phase difference Δe is reconstructed from the flow Fe by inverting Eq. (120), which yields the two following possible solutions:
We must choose one of the solutions, and we keep track of this choice by decomposing the edge set as E=E+∪E−, where E±={e∈E|Δe=Δ±e}. Not all solutions obtained in this way are physically correct. We can get the physically correct ones by making sure that the sum of the phase differences around any fundamental cycle yields 0 or an integer multiple of 2π, which is expressed via the following winding number condition:
Fixed points with E−=∅ satisfy Eq. (114) and are always stable, while states with E−≠∅ are typically unstable (Manik et al., 2014; Delabays, Coletta, and Jacquod, 2017; Manik, Timme, and Witthaut, 2017).
The cycle flow approach yields numerous analytical insights into the occurrence of multiplicity of sync states:
• No multiplicity in trees.—In a lossless tree network, either there is no fixed point or there are 2N−1 fixed points, of which one is stable and 2N−1 are unstable (Manik, Timme, and Witthaut, 2017).
• Similarly, multiplicity is ruled out in dense networks, as the fundamental cycles are too small (Taylor, 2012). Hence, all fixed points except one include edges with |δ∘i−δ∘j|>π/2.
• In a simple ring network, one can derive explicit upper and lower bounds for the number of fixed points. Generally, the number of fixed points increases with the size of the ring and decreases with the loading of the grid (Ochab and Gora, 2010; Delabays, Coletta, and Jacquod, 2016; Manik, Timme, and Witthaut, 2017).
• In plane networks, one can show that the winding numbers are unique. That is, two fixed points cannot have the same winding numbers and the same decomposition E=E+∪E−. This can be used to derive upper bounds for the number of stable fixed points (Delabays, Coletta, and Jacquod, 2017; Manik, Timme, and Witthaut, 2017).
• One can derive estimates for the number of stable fixed points using scaling relations (Manik, Timme, and Witthaut, 2017); an example is shown in Fig. 13. Further numerical results for large sparse networks were discussed by Mehta et al. (2015) and Xi, Dubbeldam, and Lin (2017).
Multiplicity of sync states: scaling of the number of stable fixed points in a grid with two cycles. The grid consists of two rings with n+1 nodes sharing one edge such that N=2n. Edge weights are assumed to be homogeneous and loads vanish ( Pi=0). The analytic scaling result 0.1576(n2+2n) (solid line) matches the numerically exact results (circles) well. Adapted from Manik, Timme, and Witthaut, 2017.
The previous approach can be generalized to include Ohmic losses (Balestra et al., 2019). One again first focuses on the lines and considers the flows Fij and the losses Lij=gjk[1−cos(δ∘i−δ∘j)] as basic variables. The power balance equation is linear in these variables, and one can construct the general solution of this set of equations. Among this large set of solution candidates, the correct solutions are found by imposing the winding number conditions (124) and the additional constraint
which follows from sin2+cos2=1. One then finds that Ohmic losses can have two conflicting effects on the existence and number of stable steady states. On the one hand, high losses must be compensated for by higher power flows, which may decrease their number. On the other hand, Ohmic losses can stabilize certain solution branches and thus foster multistability.
D. Nonlinear stability and explicit synchronization criteria
The existence of a phase-locked state is a prerequisite for the stable operation of a power grid. But even in elementary networks phase-locked states and limit cycles can coexist (see Fig. 9), and whether or not the phases lock depends on the initial conditions. In the following, we review explicit criteria for synchronization based on the dynamics (111) and discuss the main factors determining a network’s synchronization capability. A comprehensive overview is provided in the review by Dörfler and Bullo (2014).
1. Sufficient criteria for dense networks
We first present a sufficient condition for dense networks, which arise naturally in power system models after Kron reduction; cf. Sec. III.C.6. We are interested mainly in second-order power grid model dynamics [Eq. (111)], but it is simpler to start with the following overdamped limit:
Dörfler and Bullo (2012b) proved that this system achieves phase cohesiveness and frequency synchronization if the coupling is strong compared to the differences in natural frequencies [ ωnati=Pi+∑jKijsin(γij)]. They quantified these two opposing forces in terms of the parameters
Furthermore, they defined ζmin∈[0,π/2−γmax) and ζmax∈(π/2,π] as the solutions of the equation sin(ζmin)=sin(ζmax)=cos(γmax)Γcrit/Γmin whenever this solution exists.
Then they derived the following theorem: If Γmin>Γcrit and the angles are initially not too different [all δi(0) are in the arc A(ζmax)], the network will achieve exponential frequency synchronization; i.e., all frequencies ˙δi converge exponentially fast to a common frequency Ω. The phases remain cohesive and all δi(t) reach the arc ¯A(ζmin).
Dörfler and Bullo further generalized this result to the second-order power system model (111) if the damping was high enough. More precisely, if the ratio ε=Jmax/Dmin is lower than a critical value ε∗, the solution of the second-order model is well approximated by the solution of the first-order Kuramoto model up to an error of the order of ε. They further argued that, for actual power systems, ε is of the order of 0.1 if one takes into account the control system in the effective damping constant Di.
We now sketch the main idea of the proof of these theorems. For each point in time t we identify the nodes m and n with the largest and smallest angles. That is, all angles δi(t) are in the arc [δn(t),δm(t)], and
denotes the length of the arc. One can then compute the derivative dV/dt explicitly and finds that this expression is strictly nonpositive if
If this condition is satisfied, the arc A(ζ) cannot grow and the angles remain cohesive. Frequency synchronization is shown in a second step by differentiating Eq. (126) with respect to time as follows:
Under the condition of phase cohesiveness, the matrix ˜Λ(t) is a time-varying directed Laplacian. Hence, Eq. (131) describes a contraction and all frequencies ˙δn must converge exponentially to a common value Ω. Finally, the generalization to the second-order power system model is analyzed within singular perturbation theory.
2. Sufficient criteria for sparse networks
The previously presented sufficient criterion holds for all networks but is of little use if the topology is sparse since we then have Γmin=0. A different approach is needed in this case. One family of sufficient criteria for the existence of a phase-locked state uses a method introduced by Jadbabaie, Motee, and Barahona (2004). The steady-state condition (117) is reformulated in the following vectorial form:
where ˜Λ(δ) is a state-dependent Laplacian with edge weights wij=Kijsin(δi−δj)/(δi−δj). One can then apply fixed point theorems to derive conditions guaranteeing that Eq. (132) has a solution. For instance, one can formulate the following condition in terms of the algebraic connectivity of the network (cf. Sec. II), which yields a sufficient condition for the existence of a synchronized state (Dörfler and Bullo, 2012a):
where λ2(Λ) is the algebraic connectivity of the ordinary weighted graph Laplacian. We note that one has to be careful about the domain when applying fixed point theorems. In particular, the sufficient conditions do not guarantee uniqueness of a solution, as discussed by Manik, Timme, and Witthaut (2017).
Another approach uses the cycle flow arguments that were introduced in Sec. III.A.6. Denoting the sine of the phase difference along the line ℓ as ψℓ, we have the general relation (36). It has been shown that the cycle flow contribution f is actually negligible in most cases such that ψ≈E⊤Λ∗P (Dörfler, Chertkov, and Bullo, 2013). Hence, a synchronous state with a phase cohesiveness ζ exists if
Equation (134) provides a rigorous sufficient condition for many elementary networks and an excellent approximate condition for actual power grid topologies (Dörfler, Chertkov, and Bullo, 2013). Recently generalized types of Lyapunov function were derived for lossless power grids (111), allowing researchers to derive sufficient conditions for the global stability of synchronous states (Schiffer, Efimov, and Ortega, 2018, 2019).
E. Probabilistic analysis of nonlinear stability
The synchronous state of a power grid can be globally stable in favorable cases, but typically several different asymptotic attractors exist. This includes limit cycles, as shown in Fig. 9 for a single generator. In this section we explore the nonlinear stability of different attractors of the swing equation using the probabilistic methods introduced in Sec. V.C.3. We note that the swing equation is a simplified model of power system dynamics in the vicinity of the synchronized state. Hence, nonsynchronous attractors of this model do not faithfully describe asymptotic states of real systems. However, the existence of these states shows that those aspects of the system well described by the swing equation are not sufficient to guarantee stability. This includes in particular inertia and primary droop control.
1. Limit cycles and other attractors
A large variety of only partially understood attractors exist in power system models. Limit cycles similar to the ones discussed for a single machine in Sec. V.C.2 are among the most common nonsynchronous attractors. These states can be understood as arising from a partial decoupling limit. Consider a network consisting of two disconnected parts Vl and Vr. The two parts can synchronize independently at respective frequencies Ωl and Ωr according to Eq. (116).
If the two parts are weakly coupled, the equations of motion include terms proportional to sin[(Ωl−Ωr)t], which average out in time. The limit cycle with separate frequencies Ωl and Ωl persists, but the coupling causes perturbations leading to oscillations around the mean frequencies (Gelbrecht, Kurths, and Hellmann, 2020). When the coupling increases, the basin size of these attractors (the variance of the frequency time series) may even increase until the attractors eventually become unstable. Figure 14 illustrates the nonlinear stability of such states in a network with two populations of oscillators with power injections Pi=±1, effective damping Di=0.1, and homogeneous line parameters K. For weak coupling the oscillators always rotate freely at ˙δ=Pi/Di. When the coupling increases, a large variety of attractors come into being, most of them with negligible basin size (referred to as outliers). Two partially synchronous states of the type discussed earlier occur more often. The basin size first increases with K, but the fully synchronized state then becomes dominant. Its basin size increases rapidly for K≥4, reaching unity for K≳8.
Nonlinear stability of different attractors in a small sample network. (a) Structure of the network with two classes of nodes with injections Pi=±1. (b) The basin bifurcation diagram for the system with the synchronized, partially synchronized, and synchronized regimes. From Gelbrecht, Kurths, and Hellmann, 2020.
The occurrence of similar states was already studied from the point of view of hysteresis by Olmi et al. (2014). However, not all asymptotic states with an appreciable basin have this form. Olmi (2015) found that complex chimeras are possible, and Nitzbon et al. (2017) saw that perturbations at some types of nodes lead to limit cycles with a frequency not given by Eq. (116). The impact of noise, inertia, and network topology on different limit cycles and the hysteresis behavior was studied by Tumash, Olmi, and Schöll (2018).
2. The impact of network structure
We now return to the synchronized state and investigate its nonlinear stability with respect to localized perturbations in terms of the single-node basin stability. In particular, the set of initial states A is chosen such that all nodes except one are at the fixed point values initially. This measure is of high practical relevance, as actual failures or disturbances are typically limited to a single machine or a small subnetwork. Further, it reveals the vulnerable nodes in a grid. A numerical analysis for sparse networks (Menck et al., 2014) showed that the local network features have a substantial impact on the single-node basin stability. In particular, dead ends tend to cause exceptionally poor single-node basin stability. Nodes at which the dead end is connected to the rest of the grid must be viewed as the most vulnerable spots with respect to dynamical perturbations.
The more detailed study of Schultz, Heitzig, and Kurths (2014a), based on a more realistic network ensemble introduced by Schultz, Heitzig, and Kurths (2014b), identified detours as stabilizing features and showed that there is a statistically robust relationship between network features and node robustness. Kim, Lee, and Holme (2015) showed that the community structure of the system determines which nodes reach a high single-node basin stability first when the coupling strength is increased. Kim, Lee, and Holme (2016) conducted a systematic exploration of the impact of network motifs on single-node basin stability and found that betweenness and degree matter specifically. Kim, Lee, and Holme (2016) and Kim et al. (2018) showed that there is strong nonmonotonic behavior of basin stability with increasing coupling strength. Finally, Kim et al. (2019) introduced integrated basin stability as a way to understand the transition to global synchrony as a function of increasing coupling strength.
A more detailed picture emerges if we add survivability, the probability that a perturbation does not lead to a violation of transient operational bounds. Hellmann et al. (2016) showed a strong direct dependence of single-node survivability on the node degree, with high degree nodes being particularly unstable. Combining survivability and basin stability, Nitzbon et al. (2017) showed a detailed picture of various single-node stabilities as a function of topology. Key to this is a node classification into those in a loop and those on a tree. Further subdividing the tree nodes into leaves, inner tree nodes, and sprouts of high and low neighbor degree led to the identification of classes of nodes that show qualitatively different stability behavior; see Fig. 15. Perturbations at some of these nodes were already shown to result in novel asymptotic states, mentioned in Sec. VI.E.1, that cannot be explained with a decoupling limit.
Topological properties determine the single-node basin stability and survivability. Left panel: survivability (Surv) and basin stability (BS) for nodes in a synthetic power grid network ensemble, showing pronounced differences in the behaviors of various node types. Right panel: node classification according to Nitzbon et al. (2017). Edges that are not in any cycles are trees. The nodes in trees are distinguished as roots (part of a cycle and a tree), sprouts (degree 1 nodes next to a root), leaf (degree 1 nodes not next to roots), and inner tree nodes. Sprouts are further distinguished according to their neighbor degree as dense ( kav>5) and sparse ( kav≤5). Adapted from Nitzbon et al., 2017.
3. The necessity of realistic models
When identifying the imprint of network structure on the stability properties of individual nodes, one usually proceeds by simplifying the system. To isolate the impact of topology, other factors are taken out. For example, most of the previously cited studies assumed homogeneous node parameters and neglected voltage dynamics and losses. We can use probabilistic methods to understand how much these assumptions actually alter the picture.
Wolff, Lind, and Maass (2018) showed that both inhomogeneity of the nodes and the details of the coupling of generators into the grid actually have a profound impact on the single-node stability. They evaluated the basin stability and the return times to the sync state for power grids with a fixed topology but while using a variety of models including the heterogeneity in the node and line parameters and the more detailed coupling of generators, where the generator was located at an internal node coupled by a line to the system bus. The grid was completely stable for homogeneous parameters but exhibited instability when introducing heterogeneity. The coupling of generators at internal nodes was found to stabilize the system overall; however, the location of unstable nodes also changed.
Taking higher-order internal dynamics of the generators into account for single-node stability analysis was considered by Auer et al. (2016). It was found that while single-node survivability is well captured by the swing equation (111), voltage dynamics and internal nodes can lead to additional asymptotic instabilities once outside the survivability region. The first results on basin stability of higher-order generator models were obtained by Liu et al. (2019).
Finally, most of the previously discussed works and much of the theoretical literature focus on lossless grids. Real transmission lines have a small but nonzero resistance (cf. Table I), and typical values of the parameter γ in Eq. (109) are around 0.24. While neglecting Ohmic losses is often justified for the power flow in the synchronized state, limit cycles and global stability properties are strongly affected (Hellmann et al., 2020). Figures 16(a)–16(d) show the basin of attraction of the synchronous state (peach) in the slice of phase space belonging to a single node centered on the sync state. As γ is increased from zero to realistic values, the basin first expands and then switches to the other half plane.
Main results of Hellmann et al. (2020). (a)–(d) Slices of the phase space in the dimensions belonging to a node k. Different colors show different asymptotic states. Peach is the sync state, blue is a solitary state, and green is a solitary rotating against its power infeed. (e) Average single-node basin stability for the entire ensemble of networks, for the different types of asymptotic states, as well as the number of desynchronized oscillators following a perturbation (solid purple line). Adapted from Hellmann et al., 2020.
A key role in this change is played by 1-solitary states, in which the entire network stays close to synchrony except for one desynchronized node. If the desynchronized node is rotating with the driving force ⟨˙δi⟩/Pi>0, we speak of a normal solitary. If it is rotating against the local driving force, we denote this as an exotic solitary state. Figure 16(e) shows the average asymptotic single-node basin stability for different types of asymptotic states: (i) peach for the sync state, (ii) dark blue for ordinary 1-solitary states, (iii) dark green for exotic 1-solitary states, (iv) light green and blue for more complex asymptotic states that also contain individual desynchronized nodes, and (v) gray for other asymptotic states.
The presence of exotic solitary states can be understood by considering the decoupling limit for an individual oscillator from the rest of the grid. When completely decoupled it will rotate with a frequency Ωsol=Psol/Dsol. If we now increase the coupling, that state persists but, taking losses into account, the power flow no longer averages to 0 but instead to Pe=Kesin(γe) according to Eq. (109). Hence, in contrast to the lossless case the coupling to the network will not only perturb the limit cycle but also shift its location to Ωsol=[Psol−Kesin(γ)]/Dsol. This picture was confirmed in detailed numerical experiments by Hellmann et al. (2020). Methods to restore synchronization from solitary states were discussed by Taher, Olmi, and Schöll (2019).
We stress that the aforementioned results were obtained for a particular power system model, and all models face certain limitations. For instance, the 1-solitary states are related to “pole slipping” transients in practice. They correspond to true limit cycles in the restricted model, but the dynamics is changed by the generator’s protection systems, which are not included in the model. We conclude that care should be taken when interpreting model simulation results for real-world situations. Models are useful in engineering power system stability. By construction, such models idealize certain real-world aspects and thus do not capture all details of real-world power system dynamics. In particular, many dynamical models do not account for the protection systems or changes in external conditions or control systems. One advantage of probabilistic methods for stability assessment is that they can be straightforwardly adapted to models of different complexity and scope.
VII. Structural Stability of Power Grids
Structural damages are the ultimate threat to power grid stability: “Typically, the blackout can be traced back to the outage of a single transmission (or generation) element,” according to Pourbeik, Kundur, and Taylor (2006). Such an initial outage can trigger secondary outages and, eventually, a cascade of failures bringing down a power grid entirely. One of the most important security regulations is the N−1 rule, with Wood, Wollenberg, and Sheblé, 2014 having stated that “no single outage will result in other components experiencing flow or voltage limit violations.”. But this rule is violated occasionally, making a grid vulnerable and large-scale blackouts possible.
In this section, we analyze the impact of structural damages focusing on outages of transmission elements. We mostly employ a quasistatic picture, assuming that after an outage the grid will rapidly relax to a new steady state (if it exists). But in this new state other transmission elements may be overloaded, leading to emergency shutdowns or short-circuit failures, etc.
The key questions we address in the following are (1) How do failures spread in the network? Given a transmission line failure, how and where are flows rerouted? (2) How does the entire network react to structural damages or other perturbations? Is there enough redundancy to cope with damages or do secondary failures take place? (3) Finally, how do large-scale blackouts emerge? How do cascades propagate through the grid?
A. Quasistatic analysis of line outages
1. The line outage distribution factors
We first analyze the outage of a single line within the linearized load flow approximation, or dc approximation, introduced in Sec. III.A.4. Assume that a line ℓ^=(r,s) that initially carries the flow F(0)ℓ fails. If the grid remains connected, the flow change at another transmission line e^=(m,n) is linear in F(0)ℓ,
where the factor of proportionality is referred to as a line outage distribution factor (LODF) (Grainger and Stevenson, 1994; Wood, Wollenberg, and Sheblé, 2014).
We summarize the derivation of the LODFs following Guo et al. (2009) and Strake et al. (2019). As the line ℓ=(r,s) fails, the nodal susceptance matrix changes as
where the vector νrs∈RN is +1 at position r, −1 at position s, and 0 otherwise. This causes a change of the nodal phase θ→θ′=θ+α. Subtracting the linearized load flow equations (32) for the perturbed and unperturbed grid then yields
which can be solved for α and used to compute the change of the line flows as
with the superscript * denoting the Moore-Penrose pseudoinverse. Equation (137) can be greatly simplified using the Woodbury matrix identity, which finally yields
and thus LODFe,ℓ=bmnν⊤mnB∗νrs(1−brsνrs⊤B∗νrs)−1.
This procedure is readily generalized to multiple line outages (Güler, Gross, and Liu, 2007; Kaiser, Strake, and Witthaut, 2020). Assume that the lines ℓ1,ℓ2,…,ℓM fail, but the network remains connected. We define a projection matrix from the space of all links onto the subset of failing links P∈RM×L as follows:
We define projections of the node-edge incidence matrix, the branch reactance matrix, and initial flow vectors as
One can then proceed as for a single failing line. Using the Woodbury matrix identity, one obtains
In Eq. (139), we only have to invert an M×M matrix in addition to the inversion of the initial matrix B∗. If M is small, the inverse can be evaluated explicitly and we obtain a set of generalized LODFs.
We note that contingency analysis via LODFs can be improved by a modification of the linearization procedure. One starts with the nonlinear expression Fj→k=bjksin(θj−θk) for the real power flow and carries out the linearization only at a later stage (Jung and Kettemann, 2016; Manik et al., 2017). One then obtains the same results as given previously, but with reweighted line susceptances
Related approaches are sometimes referred to as a “hot-start dc approximation.”
2. Spreading of failures
A deeper physical insight into the network flow rerouting problem is obtained by an analogy to discrete electrostatics. Using again the Woodbury matrix identity, Eq. (136) can be recast in the form
with the source term
As the matrix B is a Laplacian matrix and the right-hand side is nonzero only at positions r and s with opposite sign, Eq. (141) is a discrete Poisson equation with a dipole source and α is a dipole potential; see Biggs (1997). Indeed, if we compute the impact of the outage of a line in a rectangular square grid, we simply recover a discretized version of the familiar electric dipole field, as illustrated in Fig. 17. Thus, LODFe,ℓ depends on the orientation of the two lines e and ℓ and decreases as (distance)−2.
Impact of a link failure in a square lattice with uniform edge weights. (a) Normalized change of the nodal potentials αn for a single failing link located in the center of the network. Both the size of the nodes and the color code represent αn. (b) Normalized change of the link flows ΔFe for the same topology. Arrows and color represent the direction and strength of flow changes, respectively. From Strake et al., 2019.
Understanding the impact of line outages for arbitrary grid topologies is much more challenging. The main complexity arises from the network topology encoded in the Laplacian B, which can be highly irregular. Still we can get some insight using tools from algebraic graph theory (Strake et al., 2019; Kaiser and Witthaut, 2021a).
Before we proceed to spatial failure spreading, we have a more detailed look at the right-hand side of the Poisson equation (141). The strength of the dipole |q|, or more precisely the factor χrs, can be understood from the grid topology: If a unit power is injected at node r and withdrawn at node s, then χrs is the flow via the direct link (r,s) and 1−χrs is the flow via other nondirect pathways. Hence, we can view 1−χrs as a measure of redundancy of the line (r,s). Indeed, it can be rigorously related to topological redundancy measures (Strake et al., 2019; Guo et al., 2021). An alternative approach to line outage problems was introduced by Ronellenfitsch et al. (2017). Flow changes due to a failure can be decomposed into cycle flows as ΔF=Cf using the edge-cycle incidence matrix (5). The cycle flow strength f is again determined by a discrete Poisson equation, now formulated on the dual graph.
3. Locality and the importance of distance
Intuitively, a transmission line outage should affect nearby lines more heavily than remote ones. In a regular square lattice, flow changes decay as (distance)−2, as previously discussed. But how can we quantify this statement and what does nearby mean in a network with complex topology? These questions were studied using a spectral approach by Labavic et al. (2014), Jung and Kettemann (2016), Kettemann (2016), and Rohden et al. (2016).
The discrete Poisson equation (141) can be formally solved by decomposing the Laplacian B into its eigenvectors defined as BΦj=λjΦj with indices j=1,…,N. Assuming the failure of a line (r,s), we obtain
using the generalized hot-start linear response with edge weights (140); see also Haehne et al. (2019). Similarly, the flow changes at the remaining links are given by
For a regular grid of degree 2d, the eigenvectors Φj are plain waves and the eigenvalues λj can be obtained exactly, yielding (Kettemann, 2016)
In Eq. (145) cd are constants. For d=2, a regular square grid, the power law decays with exponent 2 and c2=2π. For a random graph one finds exponentially localized eigenstates, and thus an exponentially decaying response (Torres-Sánchez, Freitas de Abreu, and Kettemann, 2020), so the disturbance due to the outage affects the power grid in only a finite range around the outage (Kettemann, 2016). This phenomenon is well known and originates from random scattering from nodes with random degree (Hata and Nakao, 2017) and/or random parameters (García-Mata et al., 2017), so-called Anderson localization (Anderson, 1958).
Numerical calculations based on the ac load flow equations were found to be in good agreement with Eq. (143); see Jung and Kettemann (2016). That study considered the addition of a new transmission line and confirmed the generic power-law decay of the change in transmitted power ΔF∼(distance)−2 up to finite size corrections. As a more realistic grid topology, the effect of a change in transmission line capacity was studied in a model of the German transmission grid (Medjroubi, Matke, and Kleinhans, 2015) that is shown in Fig. 18 (left image). The largest 2-connected component of this grid was considered, with N=260 nodes and L=479 edges, homogeneous line parameter bij and power injections randomly chosen from a binary distribution, Pi∈{−P,P}. Distance is defined as the geodesic distance; see Sec. II. The resulting response behavior depends strongly on the location of the added line. Most locations result in a long-range power-law decay with an exponent that is, on average, close to 2. However, there are certain classes of locations, denoted in Fig. 18 as subsets 1 and 2, where the change results in an exponential decay beyond a certain length scale. Notably, the regions around subset 2 are weakly connected to the remaining grid, forming weakly connected islands, which explains the fact that the disturbance decays exponentially at distances beyond these islands.
Left panel: German transmission grid topology (220 and 380 kV) (Medjroubi, Matke, and Kleinhans, 2015). Right panel: Double-logarithmic plot of transmitted power change ⟨|ΔF(ℓ)mn|⟩(r) as a function of distance r when a transmission line (ℓ) between 880 different pairs of nodes i and j was added. We averaged this quantity over all edges with the same distance r to the added line (ℓ) and performed an ensemble average over R=100 realizations with randomly distributed generators and loads. The power-law expected for a two-dimensional regular grid ΔF∼r−2 is shown for comparison (black line). The colored data belong to special subsets of added edges, as discussed in the text. Error bars, 95% confidence level. Adapted from Jung and Kettemann, 2016.
A different approach to understanding the spatial pattern of flow rerouting is to adapt the measure of distance to the respective problem (Strake et al., 2019; Kaiser and Witthaut, 2021a). Consider the line outage distribution factors LODFe,ℓ for the IEEE 30-node test grid shown in Fig. 19, fixing a failing link ℓ in the right part of the network. The LODFs decay with the geodesic distance of two links, but the correlation is only moderately pronounced. This is not surprising, as the geodesic distance does not take into account the fact that flow rerouting requires two paths between the links ℓ and e, one to and one fro. We can define an alternative distance measure, the rerouting distance, that seeks the shortest path from one end of the link ℓ via e to the other end of the link ℓ. Numerical studies show that the correlation of rerouting distance to the LODFs is much stronger than for the geodesic distance (Strake et al., 2019).
Localization of the impact of line failures as seen with different distance measures. (a) Flow changes after the outage of a single line (in red) in the IEEE 30-bus test grid. (b) Geodesic distance to the failing link. (c) Rerouting distance to the failing link accurately predicts the spatial distribution of the flow changes. (d),(e) Flow changes |ΔFe| vs distance to the failing link ℓ on a logarithmic scale. One clearly observes that the rerouting distance used in (e) yields a higher correlation than the ordinary geodesic distance shown in (d). Adapted from Strake et al., 2019.
Finding a perfect measure of distance requires one to take all details of the Laplacian B into account, not just the single shortest path. Indeed, the resistance distance (Klein and Randić, 1993) provides such a measure of distance and can be applied to line outages (Tyloo, Pagnier, and Jacquod, 2019). However, as this distance is calculated from the full matrix B, it does not offer a direct topological interpretation as the shortest path distance does.
4. The importance of community structures
Some power grids show a pronounced community structure. In the Scandinavian grid, for example, Finland is connected to its neighbors via only two ac lines [Fig. 20(a)]. This weak topological connectivity also induces a weak electric connectivity: LODFs between two lines in different communities are significantly smaller than expected from their distance. This aspect is exemplarily illustrated for the failure of a line in Sweden [Fig. 20(b)]. For most lines in Finland the respective LODFs are of the order of 10−7, but some lines in the north near the connection to Sweden are slightly more affected. A more quantitative analysis is provided in Fig. 20(c) showing that LODFs between different communities are significantly suppressed, up to several orders of magnitude, compared to LODFs within a community. Hence, a splitting into communities can strongly impede the spreading of failures (Kaiser et al., 2021). Several other topological structures have similar or even stronger effects on failure spreading (Kaiser, Latora, and Witthaut, 2021; Kaiser and Witthaut, 2021a): If the coupling between two parts of a grid can be described using an adjacency matrix of unit rank, then all mutual LODFs vanish exactly and the spreading of failures and cascades is suppressed entirely. The simplest realization of such a unit-rank coupling is given by a bridge (Ronellenfitsch et al., 2017; Guo et al., 2021) that is exploited in recent proposals for emergency measures to contain cascading failures (Zocca et al., 2021; Bialek and Vahidinasab, 2022). These results emphasize the importance of structural features for the general robustness of electric power grids.
The importance of community structures in the spreading of failures. (a) Community structure of the Scandinavian power grid obtained by spectral clustering. (b) Magnitude of the LODFs for a failing link at the border between southern Norway and Sweden (red line). (c) Magnitude of the line outage distribution factors |LODFjℓ| as a function of the unweighted edge distance of the transmission lines j,ℓ. The cases in which j and ℓ are in the same or different communities are compared. The line gives the median, while the shaded area gives the 25%–75% quantile. The grid data were taken from Hörsch, Hofmann et al. (2018).
5. Multiple outages and collective effects
When two lines k and ℓ fail, we cannot simply superpose the flow changes but instead have to account for collective effects. To see this, assume that k and ℓ fail successively. In calculating the impact of the failure ℓ, we must take into account the observation that the flow Fℓ has already changed due to the failure of k. The correct result for the impact of a multiple line outage was derived in Sec. VII.A.1 and, for two lines, reduces to
using the shorthand Leℓ=LODFe,ℓ and denoting the failing lines in the superscript. Collective effects can have surprising consequences. The first impact of a line failure may be partly compensated for by the intentional removal of a second line. Thus, suppose that the failure of line ℓ has led to an increase of the flow on a critical line e. In many cases we can find another line k such that
This concept was previously discussed by Motter (2004) and Witthaut and Timme (2015) for different types of flow networks; we return to this in Sec. VII.D. More surprisingly, we can find cases where the collective impact of two outages is the opposite of the impact of two isolated outages, i.e.,
An example of this effect is shown in Fig. 21, and a more detailed discussion was given by Kaiser, Strake, and Witthaut (2020). Further analysis shows that collective effects are particularly strong whenever the mutual LODFs, or more precisely the expression √LkℓLℓk, are large.
Collective effects can lead to a complete reversal of the flow changes in an N−2 failure. The color code of the lines indicates the magnitude of flow, with red indicating failing links. (a) Flows in the initial unperturbed grid. (b),(c) Flow changes after the individual failure of two links ( ℓ and k, respectively, in red). In both cases, the flow in the top left link ( e; boldface font) is greater than in the unperturbed grid: ΔF(ℓ)l≈0.3, ΔF(k)l≈0.1. (d) Flows after the simultaneous failure of both links ℓ and k. The flow in the top left link is smaller in magnitude than in the unperturbed grid and even changes its direction: ΔF(ℓ,k)e≈−0.7. From Kaiser, Strake, and Witthaut, 2020.
B. Robustness of power grids and critical infrastructures
The outage of a transmission element leads to the rerouting of currents and power flows, as analyzed in Sec. VII.A. This may eventually cause secondary failures of other transmission elements, which can result in a large-scale cascade affecting the entire grid. In this section we investigate the grid vulnerability to secondary failures and analyze the weak spots of a network.
1. Why secondary failures?
If the current on an overhead transmission line is too large, it will heat up due to Ohmic losses, which may lead to line sag and eventually a short-circuit fault, as discussed in Sec. V.B. Strict security rules are commonly implemented to avoid such faults such that the current on a line ℓ may not exceed a certain limit, the line rating |Iℓ|≤Imaxℓ. If voltages are close to the set point and losses are small, this directly translates into an upper limit for the real power flow
A violation of this limit typically leads to an emergency shutdown if the line is appropriately monitored. The line will not be damaged, but it is no longer available for power transmission. In practice, one differentiates between short-term and long-term loads of a transmission line. A higher loading may be acceptable on short terms, such that a higher line rating applies. Furthermore, not only currents but also voltages are important for grid stability. Strict security rules apply to the voltage level at every substation, which must remain within a certain interval around the grid reference level. A decreasing voltage may indicate a looming voltage instability, as discussed in Sec. V.A.
2. Critical links: A graph theoretic perspective
Consider a heavily loaded power transmission grid which is no longer N−1 secure. Some transmission lines may be critical in the sense that if they fail secondary overloads occur. But which lines are critical and which lines are prone to secondary overloads? These questions can be answered with extensive numerical simulations, but a better analytic understanding would be helpful. We expect outages of lines with a high load or flow to be more harmful, but this is certainly not the entire story. Whether the grid bears enough redundancy to cope with the failure is equally important. These aspects can be quantitatively understood using tools from graph theory (Witthaut et al., 2016).
Assume that a line (a,b) fails that initially carried the real-power flow F(0)a→b from node a to node b. After the failure, the power must be rerouted from a to b via different pathways. Graph theory now provides a necessary condition for this to be possible. The Edmonds-Karp algorithm yields the maximum flow that can be transmitted from a to b while respecting the line limits, which we call the redundant capacity Kreda→b of the line (a,b). The algorithm also yields the set of edges limiting the flow from a to b, which is referred to as the minimum cut in graph theory. Hence, we can obtain both a measure of redundancy and the locations of the bottlenecks in the grid.
In summary, a necessary condition for having a feasible power flow after the outage of line (a,b) is given by |F(0)a→b/Kreda→b|≤1. In practice, secondary overloads will typically occur earlier, as power flows are generally not maximum flow network flows but have to respect Kirchhoff’s laws. An example is shown in Fig. 22 for a semisynthetic power grid based on the structure of the Austrian transmission grid. One of the marked links carries a high real power flow of about 1.2 GW. It is uncritical because it has a high redundant capacity. The second line has a lower flow ( ≤0.9 GW) but lacks redundancy. The ratio |F(0)a→b/Kreda→b| is larger than 1 such that an outage of this line will always lead to secondary outages and a fragmentation of the grid. Applying the Edmonds-Karp algorithm directly yields the bottleneck that limits the redundancy.
(a) Real power flows Fab in MW in a test grid based on the Austrian grid topology. The two marked links are investigated in detail. (b) The ratio of the flow and the redundant capacity F(0)a→b/Kreda→b provides a powerful indicator for critical links. (c),(d) Illustration of rerouting pathways after the outage of links 1 and 2, respectively. Several rerouting pathways exist for link 1 such that the redundant capacity Kred is rather high, and a line failure is uncritical. In contrast, a failure of link 2 will cause secondary outages. The min cut for the rerouting problem consists of a single link (marked with an arrow). This bottleneck has a free capacity of 864 MW such that a rerouting is impossible. Grid data are based on the work of Hutcheon and Bialek (2013), with a slight adaptation of power injections. Austria was artificially islanded for illustrative purposes.
Extensive numerical tests for different power system models (Witthaut et al., 2016) show that the ratio |F(0)a→b/Kreda→b| is a powerful predictor for critical links. Furthermore, one can improve power system robustness by strengthening the bottlenecks (Rohden et al., 2017).
3. Generation variability and critical fluctuations
Thus far we have investigated potential overloads and cascades assuming that the power injections are fixed and perfectly known. But even if we know the scheduled values, there are always some residual fluctuations of generation and demand. How do these fluctuations affect power system stability? Is it possible that the grid is safe on average but that overloads occur during stochastic peaks that may eventually trigger a cascade?
These important questions were addressed by Nesti, Zocca, and Zwart (2018) using a quasistatic picture. The power injections P are treated as Gaussian random variables with a mean μP and a correlation matrix εΣP. The essential quantities in a security assessment are the relative line loads Lℓ≔Fℓ/Fmaxℓ, which are determined by the power injections via
using Eq. (33). Hence, the line loads L are also Gaussian random variables with mean μL=ΥμP and correlation matrix εΣL=εΥΣPΥ⊤. A line ℓ is overloaded if Lℓ>1. In the limit of small fluctuations ε→0 one can derive the following probability for this event using large-deviation theory:
where σ2ℓ=(ΣL)ℓℓ. While this is an approximation for nonvanishing fluctuations, it was found useful for ranking purposes. The smaller the ratio (1−|μL,ℓ|)2/σ2ℓ, the more likely an overload of a line ℓ is.
An example of such a probabilistic contingency analysis is shown in Fig. 23 for the SciGRID model of the German power grid. The nominal power injections were determined via an OPF, and the correlation matrix ΣP was derived from a statistic analysis of actual generation time series. One finds that the vulnerable lines where P(|Lℓ|≥1) is largest do not coincide with the most heavily loaded lines. High loads are one important indicator for vulnerability, but the network structure and the generation variability enter the overload probability on an equal footing via the effective variance σ2ℓ.
Fluctuations of power generation and consumption can induce overloads of transmission lines. Left panel: nominal line loads |μL,ℓ|=|⟨Lℓ⟩| in the SciGRID model of the German power grid. The nominal operation is determined from one-hour averages of the renewable feed-in and load. The feed-in of the dispatchable power plants and the power flows are then computed via a linear OPF; see Sec. III.B. Right panel: probability of a transmission line overload P(|Lℓ|≥1) when the true power injections fluctuate around the nominal values. Vulnerable lines are identified using small values of the indicator (1−|μL,ℓ|)2/σ2ℓ. From Nesti, Zocca, and Zwart, 2018.
One can extend this analysis and ask which power injections typically lead to the failure of a single line. In mathematical terms: What is the conditional probability distribution of P given that |Lℓ|≥1, i.e., given that the line ℓ is overloaded? For ε→0 this probability distribution gets sharply concentrated around the vector
where sgn denotes the sign function and eℓ is the ℓth unit vector. That is, there is a typical injection pattern leading to an overload of a line given that fluctuations are small but nonzero (which is a plausible assumption for large power transmission grids). An important conclusion is that a line failure is not necessarily triggered by large fluctuations in the neighborhood but rather can be triggered by the cumulative effects of small unusual fluctuations in the entire network. Correlations of the power injections, such as by a sudden large-scale increase in wind power generation, play a key role in such an outage.
C. Cascades of failures and large-scale blackouts
The outage of a transmission line leads to a rerouting of power flows that may cause secondary overloads, as previously discussed. In this section we discuss how this mechanism leads to a cascade of failures and how such a cascade propagates through the grid.
1. Simulating cascading failures
Cascading failures are often simulated in a quasistatic picture. One computes a series of stationary states of the grid by solving the ac [Eqs. (25)] or dc [Eqs. (32)] load flow equations, where the topology changes as elements of the grid fail. The basic steps of a cascade may be simulated using Algorithm 1.
More advanced simulations can include aspects of voltage stability as well as actions by the grid operators such as load shedding, the immediate disconnection of consumers to reduce the grid load.
We note that a wide body of literature on cascading failures exists in the statistical physics literature. Purely topological approaches based on percolation theory are particularly popular, as they allow for analytic solutions; see Albert, Jeong, and Barabási (2000), Albert, Albert, and Nakarado (2004), and Newman (2012). However, the applicability of topological models to power grid stability is limited, as discussed by Hines, Cotilla-Sanchez, and Blumsack (2010), Bompard, Luo, and Pons (2015), and Korkali et al. (2017). Topological models and measures can provide some general insight into the vulnerability of networks [see Galindo-González, Angulo-García, and Osorio (2020)], but they may be misleading when applied to particular grids. In particular, they capture neither the physics of power flows nor the heterogeneity and localization of power generation and consumption.
2. Local versus nonlocal propagation
How do cascading failures propagate through a power grid? Do subsequent failures occur in close proximity, or can they jump to remote areas of a grid?
A prime example of a locally propagating cascade was observed during the 2006 western European blackout (Union for the Coordination of Transmission of Electricity, 2007). The cascade was triggered by the shutdown of a line and a switching event in northwestern Germany. In every step of the cascade, flow was rerouted to the next available route in the southeastern direction, thereby causing secondary failures along these routes.
However, secondary outages can also take place at rather long distances from the initial failure. A simple reason is that small flow changes can also cause a line outage if a line were already heavily loaded before. In particular, a failure of line ℓ causes a secondary overload of line e if (assuming that F(0)e>0 without loss of generality)
Hence, a secondary outage occurs if either the flow change ΔFe=LODFe,ℓF(0)ℓ is large or the line e was heavily loaded before the outage such that Fmaxe−F(0)e is small. While the flow changes ΔFe are strongest in the vicinity of the failing line ℓ, little can be said about the initial line loadings F(0)e. They are determined by the grid topology as well as the power injections P and thus change every hour. Hence, the weak spots of a grid also vary.
The situation becomes more involved if we go beyond the linear dc approximation. For instance, voltage instabilities played an important role during the blackout in the western USA in July 1996, where strongly nonlocal impacts were observed (Venkatasubramanian and Li, 2004; Hines, Dobson, and Rezaei, 2017).
3. Influence graphs
Influence graphs were proposed by Hines, Dobson, and Rezaei (2017) to describe the propagation of cascades and to understand nonlocal effects. This description was based on the statistical analysis of a large number of cascading failures simulated with a detailed physical model of the grid. The simulation data are used to determine the conditional probability Rj,i,m that a line j fails in generation m+1 of a cascade given that line i failed in generation m. It turns out that probabilities do not change much with m as long as m≥1. Only the first generation m=0 is different, as the grid in the initial state is typically N−1 secure such that failures are less probable. Hence, cascading failures are readily characterized using two matrices R0,R1∈RL×L with elements Rj,i,0 and Rj,i,1+, respectively. An elementary example of a power grid and the reconstructed influence graph is shown in Fig. 24. A refined treatment including multiple line outages as separate states in a Markov chain model was proposed by Zhou et al. (2020).
Influence graphs are used to describe the propagation of cascading failures and to identify critical infrastructures. Left panel: power flows in an elementary test grid consisting of six nodes. Right panel: influence graph summarizing the likelihood that the failure of one line will induce the failure of another line. The nodes of this influence graph correspond to the lines in the original grid as identified by their end points. Adapted from Hines, Dobson, and Rezaei, 2017.
The influence graph can be used to identify critical lines as follows. If Pi,m is the probability that line i fails in generation m, then in the following generation
Upon the introduction of a vectorial notation Pm=(P1,m,…,PL,m)⊤, the probabilities evolve as P1=R0P0 and Pm+1=R1Pm for all the following steps m=1,2,…. The probability that a line j will fail at some time during a long cascade is thus given by
We can now analyze which elements have the largest impact on cascading failures. To this end one can calculate how the cumulative failure rate Ptot=∑jP∞,j will change if a single row i of the matrices is R0 and R1 is modified, representing an upgrade of the respective line i. One finds that most upgrades have little effect, but some lines lead to a drastic change of Ptot. The corresponding lines represent the critical lines of the grid.
4. Statistical analysis of cascading failures
Sections VII.C.1–VII.C.3 have provided a microscopic picture of how flows are rerouted and how cascades propagate in electric power grids. But how do these aspects manifest in real-world large-scale grids under various operating conditions? A large-scale statistical analysis of these questions was recently provided by Yang and Motter (2017) using a high quality model of the three North American synchronous power grids (Eastern, Western, and Texas Interconnection). The utilized cascade model was significantly extended in comparison to the elementary model of Sec. VII.C.1, including load shedding, generator adaptation, and a model for transmission line overheating. A large number of cascades were simulated for a variety of different snapshots of the grid: different points in time with generation and load patterns. All cascades were triggered by the simultaneous failure of three lines randomly selected from the entire grid (Western, Texas) or from a certain grid area (Eastern Interconnection). Yang and Motter (2017) then evaluated the probability P(p)ℓ that a line ℓ will fail during a cascade as well as the probability P(s)ℓ that a line ℓ will be disconnected from the grid while carrying no load after the cascade.
The first important result of the statistical analysis is that coreness is an important factor that determines the probabilities P(p,s)ℓ; see Fig. 25. Links with a higher coreness are more strongly connected and thus offer more options for flow rerouting. As a consequence they are more likely to suffer failures (the fraction of links with P(p)ℓ>0 increases with the coreness) and they also have a higher average value of ⟨P(p)ℓ⟩. In contrast, the probability of becoming disconnected (i.e., the fraction of links with P(s)ℓ>0) decreases with the coreness. Links with coreness 1 can be disconnected by a single failure, while links with higher coreness have multiple connections to the grid and are thus hardly disconnected. Notably the average ⟨P(s)ℓ⟩ still increases with the coreness. In this analysis it should be noted that the vast majority of all links have a coreness of 2.
Statistical analysis of cascading failures. The coreness of a link determines its vulnerability to fail (upper panels) or to be disconnected (lower panels) in a cascade of failures. Shown is the fraction of lines that (c) fail or (e) get disconnected in at least one cascade. (d) Probability of failure ⟨P(p)ℓ⟩. (f) Probability of disconnection ⟨P(c)ℓ⟩ averaged over all cascades and operating conditions. Simulations were carried out for the three North American grids under various operating conditions; see the text for details. Adapted from Yang and Motter, 2017.
In a second step, Yang and Motter (2017) investigated the vulnerable lines of the grids, i.e., lines with a high probability of failure ( P(p)ℓ>0.0005). It was found that the set of vulnerable lines is surprisingly small, including, for instance, only 48 out of 7637 lines in the Texas Interconnection. Moreover, there is a significant overlap of the vulnerable sets between different grid snapshots. That is, generation and load patterns determine the vulnerability only to a limited extent: many lines are vulnerable for various operating conditions. These findings raise the hope that grid stability can be significantly extended with comparably few grid extensions.
Finally, we are led to the question of which trigger events are particularly critical to grid stability. To answer this question Yang and Motter (2017) analyzed all cascades causing a load shedding above 300 MW. It was found that in these cases the three trigger events are typically close to each other and close to the vulnerable set of the grid, in terms of both geodesic network distances and geographic distances. This corresponds to our prior findings that flow rerouting is predominantly local.
5. The size of cascading failures
Cascades of failures can cause large-scale power outages. Statistics of the size of power outages are available for a variety of power grids, as shown in Fig. 26 for the North American power system. One observes that large-scale outages are actually not rare events; that is, the likelihood P of an outage vanishes slowly with the size of the outage Nout. Indeed, it has been claimed that the distribution shows a power-law behavior
for large Nout with an exponent 1.3≤β≤2.0 for various grids (Dobson et al., 2007). While it is generally hard to strictly establish a power law from empirical data (Clauset, Shalizi, and Newman, 2009), the statistics indicate that large-scale outages occur regularly.
Statistics of power outages in North American power grids (1984–1998). The empirical probability of a power outage as a function of the outage size on a logarithmic scale. The data suggest a power-law distribution with exponent 1.3≤β≤2.0. From Dobson et al., 2007.
A detailed analysis of cascading failures including empirical statistics was provided by Dobson et al. (2007). They analyzed outage statistics both in detailed numerical models and in a coarse-scaled model admitting a closed form solution (Dobson, Carreras, and Newman, 2005). This model reduces a cascade to the following three essential elements:
• Many elements with heterogeneous loading.—The model considers N infrastructure elements whose initial load L0 is drawn at random from a uniform distribution in the interval [0,Lmax]. For simplicity we normalize all quantities such that Lmax=1.
• Cascade mechanism.—If the loading of an element exceeds a limit L≥Lfail, then this element fails and the loading is distributed to other elements. For simplicity we assume that the loading of all other elements then increases by the value ΔL1.
• Initial trigger.—The cascade is triggered by some initial failures, thus increasing the initial load of all elements by an amount ΔL0.
For Lfail=1, ΔL0+ΔL1N≤1, the number of failing elements Nout follows a quasibinomial distribution
If ΔL1≪1/N, this distribution strongly peaks around the mean. But if ΔL1 increases and approaches the value 1/N, the system approaches a critical point with a power-law distribution
thereby providing a fair fit to the empirical results.
But why should a power system remain critical over many years? Dobson et al. (2007) argued that criticality emerges in a self-organized way. The loads in the system are increasing year to year; in terms of the model Lfail decreases and ΔL1 increases. As the system approaches the critical point, the probability and size of power outages increases and countermeasures are implemented, such as via an extension of infrastructures. Hence, the system remains close to the critical point. This behavior was reproduced in a detailed numerical model (Carreras et al., 2002) in which power dispatch and cascades occur on short timescales (cf. Secs. III.B and VII.C.1) and load and infrastructures evolve on much slower timescales. In the model, the slow dynamics brings the power system close to criticality and such large-scale cascades become likely. Note that the hypothesis of self-organized criticality is highly controversial, as summarized by Fairley (2004). Different explanations were put forward by Nesti, Sloothaak, and Zwart (2020) relating power-law distributions in outages sizes to power-law distributions in city sizes, and by Kaiser and Witthaut (2021b).
6. Transient effects in cascading failures
Our previous analysis of cascading failures was restricted to a quasistatic picture. We considered the existence and the properties of the steady state after the failure of a transmission line while neglecting any dynamic. However, the existence of a stable state after a line failure is only a necessary condition for grid stability, not a sufficient one. It is not a priori clear that the grid will end up in this state, but there are also other possibilities: (i) It is not guaranteed that the system relaxed to a fixed point after failure. Instead, it can also relax to a limit cycle or another attractor, as discussed in Sec. V.C.3. (ii) Even if the system relaxes to the desired steady state, it may violate operational constraints during the relaxation, which might also cause secondary outages. This relates to the aspect of survivability introduced in Sec. V.C.3. As a consequence, transient effects can increase the vulnerability of grid and a quasistatic picture can miss important threats to stability (Simonsen et al., 2008).
A detailed analysis of cascading failures, including transient overloads, was provided by Schäfer, Witthaut et al. (2018). The dynamics was modeled in terms of the aggregated dynamical model of Sec. III.C.9 and a line was instantaneously removed from the simulation when the flow exceeded a line rating |Fij(t)|>Fmaxij. The impact of such a transient outage is illustrated for an elementary test grid in Fig. 27, where a cascade is triggered by the failure of a single line (2,4) [Fig. 27(a)]. In a quasistatic picture, the grid immediately jumps to a new stable steady state [Fig. 27(c)]. However, the line ratings are exceeded transiently, leading to secondary failures of line (4,5) and eventually to a cascade disconnecting the grid.
Transient effects can increase the vulnerability of a grid to cascading failures. (a) An elementary test grid with two effective generators ( P=+1.5 s−2; squares) and three effective consumers ( P=−1 s−2; circles). All lines have Kij=10/6 and Fmaxij=1. At time t=1 s, the line (2,4) fails. (b) In a full dynamical study, line ratings are exceeded transiently, leading to secondary line outages. Shown is the number of failed lines as a a function of time. (c) In a quasistatic picture, the grid immediately relaxes back to a stable steady state that respects the line ratings. Shown are the steady-state line flows |Fℓ| before and after the initial failure. (d) Time evolution of the line flows |Fℓ(t)| taking into account both the continuous dynamics and the outage of lines. From Schäfer, Witthaut et al., 2018.
Whether transient overloads are important depends crucially on grid properties, the line ratings Fmaxij, and the initial trigger that starts the cascade. The strongest impacts have been observed for heavily loaded grids, where the likelihood of secondary overloads more than doubles.
D. Braess’s paradox
The loss of a transmission line can cause a blackout, either directly, as discussed in Sec. V.A, or indirectly via a cascade of failures. But the reinforcement of a transmission line or the addition of a new line can also induce a loss of stability (Witthaut and Timme, 2012; Coletta and Jacquod, 2016). This surprising effect is demonstrated for an elementary test grid in Fig. 28. The initial grid relaxes to a synchronized fixed point, whereas synchronization becomes impossible after certain grid extensions and the synchronized fixed point has ceased to exist.
Braess’s paradox: loss of the synchronous fixed point due to grid extensions. (a)–(c) Topology of the network. The vertices generate or consume the power Pi=±P, and the transmission lines have a capacity K. We consider the cases where (b) the capacity of one line is doubles and (c) a new line is added. (d)–(f) Dynamics of the grid given by the equations of motion (111) starting from the initial state θi(0)=˙θi(0)=0. The initial grid relaxes to a synchronized fixed point, while synchronization becomes impossible after the grid extensions. The parameters are K=1.03P, D=P, and J=1. Adapted from Witthaut and Timme, 2012.
We can understand the loss of a fixed point using the cycle flow approach introduced in Sec. VI.C.3. Increasing the capacity of a line as in Fig. 28(b) will always lower the load on this specific line and thus relieve the grid. However, any fixed point must also satisfy the winding number constraint (124) for every cycle in the network. As a consequence, all flows in the grid F are affected, which may cause the disappearance of synchronized phase-locked states, a prerequisite for normal grid operation. In particular, for the example shown in Fig. 28(b) a cycle flow must be added in the counterclockwise direction to continue to satisfy the constraint (124) after the line extension. But this would increase the load on lines 4→5 and 4→8 above the upper limit K, which is not possible. The fixed point ceases to exist.
Braess’s paradox can be used to improve grid stability. In certain situations it can be beneficial to shut down a transmission line on purpose to relieve another line (Motter, 2004; Witthaut and Timme, 2015); cf. Sec. VII.A.5. Similar effects were observed for other supply networks, starting with a seminal work on traffic networks by Braess (1968).
VIII. Perturbations, Fluctuations, and Transient Dynamics
Since the ascent of renewable energy sources, accompanied by increased trading and regulatory activities, power inputs and outputs increasingly fluctuate on all timescales. In addition, consumption patterns keep changing due to the increasing use of smart devices, digitalization, and global connectivity. Such perturbations and fluctuations require analyses of the state of a power grid as a driven, nonequilibrium system where voltages, voltage angles, and flows are nonstatic and even nonstationary as they respond dynamically to time-varying signals.
In this section, we address key questions about perturbations and fluctuations in power grids. In a mathematical model for bulk grids, we illustrate in Sec. VIII.A how fluctuations of feed-in and consumed power translate to frequency fluctuations while highlighting their non-Gaussian statistics. We also highlight the impact of grid heterogeneities on non-Gaussian features and, in Sec. VIII.B, present driving response relations for networked systems that explicitly take into account any given grid topology. Beyond deriving the general form of linear response theory, we explain collectively emerging dynamic phenomena including resonances, bulk oscillations, localization, and signal propagation in Sec. VIII.C. In Sec. VIII.D, we outline how the Wentzel-Kramers-Brillouin (WKB) method helps one to predict the probability of blackouts if rare but large fluctuations kick the system out of a stable phase-locked state. For fluctuations in production and consumption occurring simultaneously all over the grid, a natural implementation in the spirit of statistical mechanics amounts to considering ensembles of power grids (Sec. VIII.E) with possible applications to obtain macroscopic average quantities such as the market volume and market costs.
A. Frequency fluctuations from time series data
Dynamic driving signals continually change the overall state of a power grid. Because of fluctuations, the collective dynamics of electric power grids thereby is not only intrinsically out of equilibrium but also intrinsically nonstationary. Fluctuations often directly modify essential parameters of input and output and make them time dependent, reflecting fluctuating power feed-in and consumption. One key example is the fluctuating power occurring as a time-dependent parameter Pi(t), such as in the second-order model (111).
To analyze how the statistics of input power fluctuations influence the grid frequency, we first focus on the bulk frequency dynamics given by Eq. (87), assuming that the power imbalance is due to the following stochastic fluctuations at the grid nodes i=1,…,N with strength qi:
In traditional analyses of power engineering, input fluctuations are either neglected or modeled to be of a Gaussian nature, with either white or colored noise, i.e., without any or with small temporal correlations; see Zhang and Li (2010), Wood, Wollenberg, and Sheblé (2014), and Schäfer et al. (2017). For Gaussian white noise power fluctuations, the resulting distribution of frequency deviations ¯ω is Gaussian with zero average and a standard deviation σ=√∑iqi/2η¯J2. Such a Gaussian model neglects heavy tails observed in bulk frequency distributions (Schäfer, Beck et al., 2018; Gorjão, Jumar et al., 2020). As Fig. 29 illustrates, such heavy tails indicate that strong deviations from the reference frequency are orders of magnitude more frequent than a Gaussian model would predict.
Non-Gaussian frequency deviations with a reduced number of extreme events. (a) Distribution of frequencies around a set frequency (50 Hz) in the Continental European grid (yellow) deviates from the best Gaussian distribution (parabolic blue solid line). In particular, the probability (per Hz) of large events of more than 0.1 Hz frequency deviations is increased by more than a factor of 10 compared to the Gaussian fit. (b) Counts of “extreme events,” i.e., the number of frequency deviations larger than 0.1 Hz resolved for each minute, within full hours, accumulated for 2017 and 2011, respectively. Compiled using data from Reseau de Transport d’Electricite (2017) and the methods of Schäfer, Timme, and Witthaut (2018).
Heavy-tailed frequency distributions may emerge from either heavy-tailed distributions of the power fluctuations ¯ξ(t) or the temporal variability of the system parameters (Schäfer, Beck et al., 2018). In particular, abrupt changes of power generation at the beginning of the trading intervals may cause large frequency deviations (Gorjão, Anvari et al., 2020). Counting the number of extreme events of the frequencies, defined as frequency deviations of more than 100 mHz from the set frequency, and comparing these data from the years 2017 and 2011 reveal additional information [Fig. 29(b)]; cf. Schäfer, Timme, and Witthaut (2018). First, the total number of extreme events is reduced in 2017. Second, the number of these events is predominantly reduced in an interval of a few minutes after full hours. Third, the total number of threshold violations seems to still be substantial in the first few minutes, as well as in about the 18th, the 33rd, and the 48th minutes after a full hour (each a few minutes after a full quarter hour), hinting that switching the trading intervals from each hour (in 2011) to each quarter of an hour (in 2017), and thereby reducing the volume of trading per event together with changes in regulatory action, might have caused the reduction in the number of extreme frequency deviations. The histogram also hints at a characteristic timescale of response of frequency deviations at approximately 2–4 min delay. Both the details of such deviations and the generality of their occurrence remain poorly understood to date.
The resulting extreme events pose theoretical questions for analysis and serious practical challenges, such as for security assessment. Recent modeling work (Wolff et al., 2019) illustrated that in power grids in which consumption, generation, and transmission infrastructures are heterogeneous, fluctuating wind power injection at nodes that are weakly coupled to the grid particularly contribute non-Gaussian features to frequency deviations. Some further studies considered non-Gaussian effects, focusing on either theoretical aspects of how they may emerge or numerical evaluations of individual wind and solar data; see Kashima, Aoyama, and Ohta (2015), Anvari et al. (2016), Schmietendorf, Peinke, and Kamps (2017), and Totz, Olmi, and Schöll (2020). In particular, the intermittent nature of short-term wind fluctuations cause novel types of frequency and voltage fluctuations and thereby influence stability properties of grid dynamics (Schmietendorf, Peinke, and Kamps, 2017). Wind power induced fluctuations, quantified in terms of frequency deviation variance, moreover, propagate along interconnected synchronous machines in a characteristic way, with exponentially decaying amplitudes (Haehne et al., 2019).
B. Network linear response theory for fluctuating power
How does the collective grid dynamics respond to input fluctuations? We review the general linear response theory valid for small perturbations with arbitrary time dependence (Haehne et al., 2019; Zhang et al., 2019; Zhang, Witthaut, and Timme, 2020). Consider a small perturbation δP(t)=[δP1(t),…,δPN(t)]T that at all nodes i results in phase deviations αi(t) To first order in the deviation αi, perturbations evolve according to the linear wave equation (112), as governed by the weighted graph Laplace matrix L∈RN×N, which is defined in Eq. (113). For example, for homogeneous parameters the inertia Ji=J, and the damping factor Di=D, the response of the phase αi(t) to dynamic perturbations, such as a change of power δPj(t) at node j, is given by (Haehne et al., 2019)
with the grid reference frequency ωR; see Sec. III.C.1. The propagator from node i to node j is defined by
See Haehne et al. (2019) for the full derivation. In Eq. (158) 1/τ=D/J is the local relaxation rate (of a single node disconnected from other nodes) and ˜Λn=(J/D2ωR)Λn, where Λn∈R are eigenvalues and ϕn∈CN are the corresponding eigenvectors of the generalized graph Laplacian matrix L [Eq. (113)], which is related to the stability matrix used in small signal stability analysis (Milano, 2010; Zhang, Rehtanz, and Pal, 2012). For each eigenvalue of the generalized graph Laplacian Λn∈R there are two eigenvalues of the linearized swing equations (84), the linear wave equation (112), as given by ϵnσ=−i(1+σ√1−˜Λn)/τ, where σ=±. Thus, the denominator in Eq. (158) is proportional to the difference of these eigenvalues 2√1−˜Λn=τi(ϵn+−ϵn−). Note that the relaxation rate of each mode is given by Γnσ=−Imϵnσ, which is identical to the local relaxation rate 1/τ if ˜Λn>1 but differs from it by the term σ√1−˜Λn/τ if ˜Λn<1, yielding for the mode σ=−1 a slower relaxation in the grid than for individual nodes. This expression applies to any disturbance δPj(t) as long as its amplitude is sufficiently small; see Eq. (S2) in the Supplemental Material of Tamrakar, Conrath, and Kettemann (2018) for a validity condition. Similar expressions can be obtained in response to a change of any of the system parameters, like a change of power capacitance δKij(t) between the nodes i and j (Kettemann, 2016; Witthaut et al., 2016; Manik et al., 2017).
In the frequency representation of the linear response, we write the phase deviation αi(t) as a generalized Fourier series and expand its spatial dependence in terms of eigenvectors ϕn of the generalized Laplacian L. Thereby we obtain (Kettemann, 2016; Auer et al., 2017; Tamrakar, Conrath, and Kettemann, 2018; Zhang et al., 2019)
where cn(ε) is the contribution strength of the angular frequency ε at the nth node. Likewise expanding the disturbance in a Fourier series, we obtain
Inserting these expansions for αi(t) and δPi(t) into the linear wave equation (112), one finds, requiring that the equation is fulfilled for each term of the Fourier series, (−τ2ε2−i2τε+˜Λn)cn(ε)=ηn(ε). For a given disturbance, the Fourier component of the phase deviation cn(ε) is thus given in response to that of the disturbance ηn(ε). Inserting that expression for cn(ε) back into the Fourier series, one gets
A linear response theory generalizing Eq. (157) to include inhomogeneous parameters and the presence of Ohmic losses was derived by Plietzsch et al. (2019), and the case without inertia ( J=0) was analyzed by Tyloo, Coletta, and Jacquod (2018). Zhang et al. (2019) disentangled the responses and analyzed spatiotemporal response patterns (to be discussed in Sec. VIII.C) starting with a focus of driving one node at one frequency, i.e., ηn(ε)=δ(ε−ε0)δn,j in Eq. (160), to obtain characteristics of linear response estimates analytically.
Given the phase deviation αi(t) for the time t and the position i, one can calculate the temporal and spatial evolution of phase deviations, the frequency deviations (δf)i(t)=∂tαi(t), and the rate of change of the frequency deviations ∂t(δfi(t)) (Pagnier and Jacquod, 2019a), as well as time-averaged moments of these quantities. Furthermore, moments of increments of frequency at the node i, fiΔt=fi(t+Δt)−fi(t) (Haehne et al., 2019), contain information about correlations at timescale Δt. In the following, we review the results obtained in linear response for these quantities for different types of nonstationary signals impinging on exemplary power grid structures.
C. Spatiotemporal responses from localized to resonant
1. Propagation of short-duration disturbances
For disturbances that last for a short time compared to the local relaxation time τ=J/D, the linear response at the node i due to a power pulse
at the time t0 and the node j by Eq. (157) yields, for the times t>t0,
For meshed, spatially embedded grids and sufficiently large inertia, analytical (Kettemann, 2016) and numerical results (Tamrakar, Conrath, and Kettemann, 2018) showed that disturbances propagate ballistically with the velocity v such that the time t>t0 when signals arrive at the geometrical distance r from the position of the disturbance is given by
In Eq. (164) the arrival time t>t0 is defined as the time that it takes after a single-node disturbance in δP at time t0 until a phase deviation αi(t) exceeds a threshold αth (in the example in Fig. 30, αth=10−6JδP/D2ωR is chosen).
Arrival time t∗=(t−t0)/τ of a short pulse disturbance in units of τ=J/D vs the geometric distance r∗=r/a in units of the transmission line length a. The green rhombi represent the numerical results for (a) a square grid and (b) German transmission grid topology: KJ/D2ωR=105, a realistic value for high voltage transmission grids. Also shown are the ballistic equation (164) with the fitted velocity v (red line) and the analytic velocity equation (165) (red dashed line). The black squares indicate node-resolved arrival times t∗i. Unfilled squares indicate arrival times, averaged over all nodes at the same distance (c) for a Cayley tree with branching b=2, and (d) for the German transmission grid. In (c) and (d) KJ/D2ωR=10. Diffusive equation (166) fitted (pink line) and with the analytical equation (167) (pink dashed line). Also shown are the ballistic equation (164) fitted (red line) and with analytic velocity (red dashed line). σ=P/Pc=0.1, with critical power Pc for the respective grids. Adapted from Tamrakar, Conrath, and Kettemann, 2018.
For homogeneous regular grids, the velocity v can be derived as follows from the response theory equation (163) (Kettemann, 2016; Tamrakar, Conrath, and Kettemann, 2018):
with the length of a single transmission line a. Here, P<Pc, where Pc is the critical power above which no stationary solution to the power flow equations exists. Equation (165) gives lower bounds for the arrival times [Eq. (164)] for both regular square and German transmission grid topology [the red dashed lines in Figs. 30(a) and 30(b)]. According to Eq. (165) the maximal velocity of disturbances in both regular and meshed real grid topologies increases with the decreasing inertia J and decrease as P→Pc.
For unmeshed grid topologies it was found by Tamrakar, Conrath, and Kettemann (2018) that arrival times grow not linearly but quadratically with distance r, resembling the diffusion
For Cayley tree graphs with branching number b the diffusion constant is derived as
where Δ=(K/JωR)1/2(1−P2/P2c)1/4(√b−1) is the positive Fiedler value, the spectral gap of the graph Laplacian, independent of the number of grid nodes N. This is confirmed by the numerical calculations illustrated in Fig. 30(c) (dashed line). For low inertia ( J<Jc) the collective dynamics of coupled nodes also results in diffusive spreading of the disturbances in meshed grids; see Fig. 30(d). Jc is obtained from the condition that slow modes with a small relaxation rate appear where the spectral gap (Fiedler value) Δ is smaller than the local relaxation rate, ( 1/τ=D/J) (Tamrakar, Conrath, and Kettemann, 2018). As Fig. 30(d) illustrates, the spreading is more strongly delayed for some nodes, and disturbances are localized in certain grid regions where nodes do not become excited above the threshold (at least within the observation time). Linear response theory may explain this feature, as the response [Eq. (163)] is proportional to the eigenvector amplitude of the Laplacian ϕni. As previously noted, the long-time transient behavior is dominated by the Fiedler vector, the eigenvector with the smallest nonzero eigenvalue, which is shown in Fig. 31 (left panel) for the German transmission grid topology and illustrates two geographic regions with high amplitudes. A similar result was obtained by Pagnier and Jacquod (2019a), who found the Fiedler vector of the weighted European transmission grid to be localized at the southern and northern borders. Moreover, the global RoCoF was found to decay with increasing system inertia, which they related to the Fiedler vector intensity. In random graphs, the Fiedler vector can be strongly localized, even on a single node [Fig. 31 (right panel)], so single-node disturbances of sufficiently small amplitudes remain localized. Thus, strong randomness and inhomogeneity may result in localization of disturbances, as noted by Kettemann (2016), a phenomenon well known as Anderson localization (Anderson, 1958).
Normalized Fiedler vector intensity |ϕ2,i|/maxn|ϕ2,n| relative to its maximum (shown in black) for the German transmission grid (left panel) featuring maximal intensity at the south and north borders and for a random grid with strong localization at a single node (black; at the bottom right of the right panel). Adapted Torres-Sánchez, Freitas de Abreu, and Kettemann, 2020.
2. Propagation of stochastic disturbances
On large timescales, frequency control measures compensate feed-in fluctuations of renewable generators, as reviewed in Sec. III.D, thereby maintaining stable grid operation. However, on timescales below 1 s, grid frequency fluctuations increase with increasing wind power production (Haehne et al., 2018). Moreover, the timescale separating local from interarea modes is also of the order of 1 s (Zhang, Rehtanz, and Pal, 2012). Are such fluctuations a local feature, for instance, resulting from locally high wind power injection, or do they affect grid dynamics over large ranges? To address this question, the subsecond grid frequency dynamics has been simulated by stochastically perturbing the grid. Model simulations of coupled nonlinear oscillator models with synthetically generated wind power feed-in time series (Haehne et al., 2019) indicate that the variance of short-term fluctuations decays for large inertia exponentially with distance to the feed-in node. These findings hold for both linear chain networks and German transmission grid topology (see Fig. 32), which is in agreement with analytical results for the variance of frequency increment distributions,
obtained by linear response theory [Eqs. (157) and (158)] for chainlike grids with N≫1 nodes. Here ^DΔtf(t)=f(t+Δt)−f(t) for any function f, and ⟨(^DΔtδP1)2⟩ is the second moment of the increment distribution of the disturbance at site j=1. Thus, the second moment of the frequency increments decays exponentially with topological distance di,j=i−1 from the position of the disturbance with correlation length
Below a critical inertia J<Jc=ωRD2N2/π2K≈1.6×106 kgm2, there are nonzero eigenvalues Λn<1; cf. Fig. 32. According to Eq. (158), the amplitude of these modes decays with a rate smaller than 1/τ.
(a) Variance of increments ⟨^DΔtω2i⟩ for a chain of N=50 oscillators. The straight lines have a slope m=−1/ξ with Eq. (169). (b) Eigenvalues of Laplacian Λn for various inertia J. (c) Fits of slope m (orange) converge to an analytical m=−1/ξ [Eq. (169)] (blue) with increasing inertia J. Vertical dashed line, Jc; error bars, 2σ confidence bound. From Haehne et al., 2019.
In sharp contrast, the kurtosis
of frequency increments, quantifying deviations from Gaussian distribution ( k=3), is found to decay slowly, subexponentially, with distance from the disturbance. Thus, the non-Gaussian shape of frequency fluctuations (see Sec. VIII.A) persists over long ranges (Fig. 33).
Kurtosis k of frequency increment PDF in a chain of N=100 nodes. The frequency increment distribution p(^DΔtωi) deforms only slowly toward an almost Gaussian distribution. Insets (from left to right): i∈{2,20,50}, θ=0.01 s. From Haehne et al., 2019.
In addition to these fundamental aspects, linear response arguments have been used in various applications. Using realistic grid models, analyses have addressed how fluctuations affect the primary control effort (Tyloo and Jacquod, 2021), where additional inertia should be placed (Pagnier and Jacquod, 2019b), and which fluctuation sources have the strongest impact on the grid (Gambuzza et al., 2017).
3. Localization, distributed resonances, and bulk oscillations
Which collective dynamical phenomena does linear response theory capture? Zhang et al. (2019) analyzed the response dependence on grid topology, on the exact location of perturbed and responding nodes in the network, and on the frequency content of the power fluctuations driving the system and compared the results to the full nonlinear system dynamics. If a given node j is driven by power fluctuations at a given frequency ω, i.e., ηn(ε)=δ(ε−ω)δn,j in Eq. (160), three regimes with qualitatively different stationary responses emerge; see Fig. 34.
Localization, distributed resonances, and bulk oscillations. In a small network example (a), the joint implications on response strength (b) of driving one node (labeled 0) are illustrated as a function of driving frequencies [(c)–(e)] and at four selected response nodes (labeled 1–4). Whereas at high driving frequencies responses are localized (algebraically in frequency, exponentially in internode distance), they are irregular across frequencies and among nodes in a resonance regime of intermediate driving frequencies. At low frequencies, globally homogeneous bulk oscillations emerge where the entire grid follows the driving signal (even if at only one node). All response strengths are quantified in terms of the amplitude of the frequency response relative to their low frequency limit as ω→0. From X. Zhang; cf. Zhang et al. (2019).
For large frequencies (compared to the range set by the Laplacian eigenvalues), the responses are strongly localized on the network, with amplitudes
in Eq. (159) asymptotically (as ω→∞) decaying exponentially with frequency and algebraically with graph theoretic distance d=d(i,j) response node i has from the driven node j. In Eq. (171) const is a constant value, independent of frequency yet containing information about eigenvector overlaps and the identity of the driving and response nodes. This result complements the findings on localization reported in Fig. 31 for short-duration perturbations to signals with arbitrary frequency content and asymptotically quantifies localization. An exact asymptotic expansion up to order ω−2d−1 shows that the entry (i,j) of all powers m<d(i,j) of the Laplacian exactly equals 0, reflecting no existing paths from node j to i of length shorter than their topological distance d(i,j); see the Supplemental Material of Zhang et al. (2019).
For intermediate frequencies in and near the range defined by the Laplacian eigenvalues, resonances occur that induce spatiotemporal responses that are strongly inhomogeneous both as a function of frequency ω and among response nodes i. Near resonance frequencies relative to response strengths substantially exceed those at low frequencies; see Fig. 34(b) as well as Fig. 2(f) of Zhang et al. (2019). Here the interaction network topology plays a major role in selecting which node and frequency responses are particularly strong (or weak), thereby stressing the distributed nature of these resonances. At low driving frequencies, the dynamics at all grid nodes follow the driving signal almost instantaneously, resulting in spatially homogeneous bulk oscillations.
Such an analysis transfers to fluctuating signals with distributed frequency content and distributed driven nodes, for instance, grid models driven by purely random processes as well as those driven with power frequency fluctuations characteristic of photovoltaic or wind power generators are well characterized by linear response theory. Only above 95% transmission line loads, a regime to be avoided in real operations for various reasons, does linear response theory yield substantial errors.
D. Blackouts as rare events due to large fluctuations
Dynamical instabilities of the power grid dynamics will become increasingly important as fluctuations from renewable resources become more frequent and possibly stronger. How close is a given system to an unstable state, and thus what is a typical timescale for fluctuations to kick the system out of its stable operating state? If strong fluctuations are relatively infrequent, it is desirable to quantify how rare they are. In addition, one would naively expect that fluctuations of power production with non-Gaussian distributions increase the risk of desynchronization, as large fluctuations are less suppressed in their broad tails. An analytic framework pursued by Hindes, Jacquod, and Schwartz (2019) quantifies the risk of escape depending on the deviations from Gaussian fluctuations. This work generalizes the analytic estimates of Schäfer et al. (2017) outlined in Sec. V.C.4 in several respects. It captures the distributed dynamics of large grids and Poissonian (rather than Gaussian) noise suited to model fluctuations in real-power production.
Hindes, Jacquod, and Schwartz (2019) employed a WKB approach applied to classical stochastic systems; cf. Dykman et al. (1994). They found that the rate of desynchronization may exponentially speed up or slow down, depending on how the statistics of fluctuations combines with the least stable mode of the network described by the Fiedler vector. In contrast to the Kramers-like formula for the escape rate [Eq. (107)], which depends only on the second moment, higher-order moments of the Poissonian noise are captured by the WKB framework, and their impact on the escape rate may be counterintuitive. As we later argue, the escape time ⟨T⟩ is proportional to exp[−(S(0)+Δ(n)S)], with the action S expanded in powers of the distance to the bifurcation point. Higher-order corrections Δ(n)S to the leading term S(0) can have both signs, increasing or decreasing the escape time. We summarize the assumptions for this WKB approach to sketch the main steps of deriving the escape time from a stable state. For detailed derivations see Hindes, Jacquod, and Schwartz (2019) and references therein.
Consider power grids described by the swing equation for synchronous machines, where an input power ¯Pi is subject to fluctuations pi(t). Given histogram data, such as from turbulence induced wind-generated power increments, we construct a time series pi(t) of fluctuations reproducing the increment histogram with bins of size b. The ansatz is overdamped dynamics
with damping rate γ, and a statistical drive ξi(t), so large intermittent spikes are allowed, as they are observed for wind and solar sources. The stochastic drive is given as ξi(t)=∑bngibδ(t−tib[n]), where the power increment gib at i is the average of the pulse amplitude over the bin b∈{1,2,…,M} of the histogram, where averaged amplitudes are the measured increments gi=pi(t+1)−pi(t) and tib[n] denotes the time at which the nth such increment in bin b at unit i occurs. For modeling, the driving signal is assumed to be Poisson shot noise, so the time between two events where unit i receives a power increment within bin b is exponentially distributed with a rate νib that is consistent with the histogram.
Consider a synchronized phase-locked state that emerges through an inverse saddle-node bifurcation as the coupling strengths increase such that unstable, saddle phase-locked states are in the vicinity of it. Moreover, take typical fluctuations to be small compared to distances to the saddles and large fluctuations to be rare. Such fluctuations are not captured in the large-deviation approach of Nesti, Zocca, and Zwart (2018) described in Sec. VII.B.3, where an overload may also be due to the many small deviations that accumulate.
The thinking behind the analytical approach of Hindes, Jacquod, and Schwartz (2019) is that a large fluctuation drives the system along a most-likely path in phase space. The path connects the stable phase-locked state with the saddle, while fluctuations and the network dynamics coact as to maximize the probability of desynchronization. To derive the optimal path from classical mechanics, a generalized Fokker-Planck equation for the stochastic network dynamics is adapted to incorporate Poissonian noise. It determines the probability distribution ρ of finding oscillator phases Φ, phase velocities v, and fluctuations p at time t. To analyze large fluctuations that are rare, we project on the exponential tails of ρ and insert a WKB ansatz
into the Fokker-Planck equation. Next the action S(X+δX) with X≡(Φ,v,p) is expanded in deviations δX about X (here the stationary state), where S(X)≫1, while |δX|∝1 is assumed and only leading-order terms in the first derivatives ∂ΦS, ∂vS, and ∂pS are retained. As a remark, this ansatz allows one to analyze the shape of the probability distribution even in the remote tails for which δX∝Ω≫Ω1/2 if Ω1/2 is the typical size of a fluctuation in the stationary state, in contrast to the typical and small fluctuations considered in Eq. (106) in Kramers’s escape theory or entering Eq. (150) as it is derived from large-deviation theory. The insertion of Eq. (173) into the Fokker-Planck equation yields a Hamilton-Jacobi equation for an “action” S(Φ,v,p,t), from which a classical Hamiltonian H(Φ,v,p,λΦ,λv,λp) can be read off. The Hamiltonian depends then on “coordinates” Φ, v, and p, and their conjugate momenta λΦ,λv,λp equates to zero [ H(Φ,v,p,λΦ,λv,λp)=0]. The action reduces to
The next goal is to determine the optimal (least action) path in phase space as a solution to the classical Hamilton’s equations of motion, subject to the boundary conditions to connect the stable fixed point with the unstable saddle. The action along this optimal path is stationary [therefore, we skipped the explicit time dependence of S in Eq. (174)] and enables us to estimate the expected waiting time ⟨T⟩ for desynchronization via
where Φs denotes the stationary solution. To estimate S(Φs,0,0), the solution is determined for coupling strengths K=KSN(1+κ) close to KSN, where saddle-node bifurcation happens. To lowest order in κ, parametrizing the distance to the bifurcation point, the solution on the one-dimensional submanifold takes the form
with x∈[−1,+1] such that, for x=−1, Φi(t)=ΦSNi−Cκ1/2ri=Φ∗i. The stable fixed point is the starting point, and for x=+1 we have Φi(t)=ΦSNi+Cκ1/2ri=Φsi as the unstable saddle as the end point of the optimal path; ri is the component i of the Fiedler mode of the Laplacian that encodes the network topology. This mode is most sensitive to external fluctuations and mediates the unlikely escape. The constant C is independent of xi and κ and depends on the grid adjacency matrix, the phases at the bifurcation point ΦSNi, and the Fiedler mode ri.
As it turns out (Hindes, Jacquod, and Schwartz, 2019), to lowest order in κ the action S(Φs,0,0) depends on the damping constant γ, the saddle-node coupling KSN, and the damping rate of power fluctuations squared (γ2). Moreover, each node of the grid i contributes to S, a term proportional to μ−12, with μ2=∑bνibg2ib the fluctuation variance. This means that to lowest order in the distance κ from the saddle, the expected time to blackout ⟨T⟩∝e−S is insensitive to higher-order moments μn, n≥3, of the fluctuations (for Poissonian noise). However, at higher orders in κ, the corrections to S and Δ(n)S scale with κ according to κn−1/2. They depend on higher-order moments μn and on the balance between the positive and negative Fiedler’s mode components. For n=3, it is the product of the skewness of the power fluctuation distribution with the skewness of the Fiedler mode that determines whether the non-Gaussianity leads to an increase or decrease of the desynchronization rates captured by Δ(n)S. Symmetric power fluctuations with μ3 close to zero yield an increase in the desynchronization rate, as naively expected. Overall, this result (Hindes, Jacquod, and Schwartz, 2019) demonstrates how the network topology interacts with the stochastic dynamics in a nontrivial way.
A further application of the relation between desynchronization and the Fiedler-mode values amounts to a dimensional reduction of the phase space. The phase space is usually high dimensional for a power grid, but the possible desynchronization paths lie in a low-dimensional subspace under the given assumptions. For networks in which the saddle-node bifurcation is induced by a single overloaded edge with phase difference π/2 [see also Fliscounakis et al. (2013), Manik et al. (2014), and Rohden et al. (2017)], a so-called synchronized subgraph approximation becomes exact (Hindes, Jacquod, and Schwartz, 2019). The Hamilton’s equations of the two subsystems that desynchronize reduce to a single noisy oscillator system in relative phase-space coordinates. More generally, the subgraph approximation is sensible if the corresponding network partition is guided by approximately uniform Fiedler-mode values. Thus, rare desynchronization events become analytically predictable in spite of the high-dimensional phase space of the grid dynamics.
E. Ensembles of power grids for distributed fluctuations
In Secs. VIII.A–VIII.C the focus was on the spatiotemporal propagation of a perturbation applied to a single or a few nodes, or the addition or removal of individual lines. Mureddu et al. (2015) introduced a suitable framework for estimating the overall energy mismatch between day-ahead estimates and real data of an entire grid that should be balanced with an appropriate trade on the energy balancing market. Differences between production and consumption result from perturbations at all nodes simultaneously. The framework was later employed in different applications (Mureddu and Damiano, 2017; Korjani et al., 2018; Mureddu and Meyer-Ortmanns, 2018).
1. Basic approach
The idea is to consider an ensemble of power grids in analogy to an ensemble of microstates in statistical mechanics. The power grids differ in their “microstates,” which enter average values of global observables. Global observables may be the market volume, the amount of energy that the energy balancing market must compensate for, the costs for this amount of energy (Mureddu and Meyer-Ortmanns, 2018), the resilience of the grid (Mureddu and Damiano, 2017), or the optimal positioning of storage devices (Korjani et al., 2018).
The individual “configurations” that represent the microstates are generated from a first reference configuration of producers and consumers that is representative of part of the considered region, a certain partition between renewable and conventional power generation, a certain time during the day, or a season of the year. The reference configuration is chosen to satisfy the OPF equations; it is regarded as the operation point of the system. For either transmission or distribution grids, the ensemble of microstates is generated by applying fluctuations in production and consumption relative to the reference configuration.
The fluctuations can be due to load fluctuations, forecast errors for renewables, intraday electricity trading, or other factors. In the simplest case, the fluctuations are chosen from a truncated Gaussian distribution. Gaussian-like forecast errors are defined by standard deviations σi at nodes i and represent the expected power variations at that node at a given time. For a load of type ℓ (where ℓ∈{w,pv,…} denotes wind, photovoltaic, and other types of power generation) and power demand Pℓ, the standard deviation of the forecasting error is denoted as σℓ, and mℓ and Mℓ are the minimum and maximum values of the support of the distribution ρℓ of Pℓ, respectively, representing power constraints of the different generators. Depending on the different assignments Pℓ to nodes i, a microstate is sampled by adding a random value to the expected power production or consumption at every node i, where the random variable is extracted from the truncated probability density function ρti(x) as follows:
where the original densities (with support in all of R) are ρi(x)=(2πσ2i)−1/2exp(−x2/2σ2i) for wind, photovoltaic, or other generation at nodes i. For skewed distributions, modeling, e.g., power production by wind, we may choose the following Weibull distribution:
with a=2 (Lun and Lam, 2000; Seguro and Lambert, 2000; Lu et al., 2013) and λ=2Pw/√π, where Pw represents the wind power production of the reference configuration and Γ denotes the gamma function. Specific values of the parameters fixing the respective distribution depend on the load and the type of renewable energy and are chosen from recorded data.
As a next step, quantities like the resulting mismatch Pji of the power at node i in configuration j are measured and the Pji summed over all nodes. This way we obtain the total mismatch Pj in power production (the so-called market volume) for a given configuration j, which the market will balance. To obtain representative values of the market volume, we sample a sufficiently large number of configurations to include in the ensemble. In particular, one can analyze how the distribution of the market volume over the volume size depends on the distribution of fluctuations (normal or Weibull), and (via the choice of the reference configuration) on the time of day, the season, the geographical zone, and the percentage of renewables with respect to the total production (Mureddu and Meyer-Ortmanns, 2018). The distribution of this market volume enters the market costs, the average prices, and the average profit per technology that is involved in the production.
2. A second network layer
To model the energy market and to determine market costs, prices, and profit, a second network layer describes market trading. Here one option is agent-based modeling with network nodes representing agents [cf. Han et al. (2019)] such that the agents are the retailers, with one retailer per conventional power station of the physical grid. In the work of Mureddu and Meyer-Ortmanns (2018), the agents first have to undergo a learning phase, in which they learn how to place optimal bids in competitive auctions with the aim of buying (or selling) in the most profitable way. To simulate how real market operators acquire knowledge about the market in time and adapt their decisions, a modified Roth-Erev algorithm (Nicolaisen, Petrov, and Tesfatsion, 2001; Mureddu et al., 2015) adjusts the offer propensities of agents in a self-consistent way, with the goal of maximizing profits. The agents interact via a so-called market authority that provides a link between the physical grid and the market. The market authority knows the mismatch from the physical grid and takes the bids of the retailers, accepts or rejects these bids, informs the retailers about the decision, and proceeds until the required mismatch in energy is covered at the lowest possible costs.
If the energy balance is restored with energy provided by the subset of retailers who offered the energy at the lowest price, the feedback on the physical grid stability must be checked since the economically best selection need not be reasonable from the viewpoint of grid stability. Lines might get overloaded if the conventional generators that were selected for selling the energy happen to spatially cluster together. Thus, this framework of different network layers, coupled via a market authority, allows one to analyze the feedback from economical (low costs) to physical (high stability) optimization objectives.
The physical grid stability is particularly endangered if the retailers behave like arbitrageurs when the reserve energy price falls below the price on the intraday market. In such a situation, retailers play a kind of minority game (Ritmeester and Meyer-Ortmanns, 2021). Modeling their behavior accordingly leads to suggestions on how to control arbitrage. For example, if the few large parties contributing to the market (rather than the many small ones) are made risk averse due to small penalties, it has a disproportionately large effect on reducing the abuse of price differences in terms of arbitrage. Moreover, from the remarkable analogy of the minority game with spin glasses it becomes understandable why a larger number of retailers may not at all reduce the fluctuations in arbitrage, as one would naively expect. Related features of underlying phase transitions are also visible in realistic markets.
3. Alternatives for treating uncertainties
A number of alternative approaches exist to deal with inherent uncertainties in power and demand. The goal is to derive distributions of induced fluctuations. Induced fluctuations may refer to the energy mismatch (as previously discussed), or to induced voltage or frequency fluctuations. An approach called chance-constrained ac optimal power flow (ac CC OPF) was described by Roald and Andersson (2018) and references therein. In contrast to the optimal power flow summarized in Sec. III.B, chance constraints ensure that the system constraints will be satisfied only with a specified probability. The framework of ac CC OPF optimizes not only the scheduled dispatch with respect to costs but also the procurement of reserve and voltage control during deviations from expected values. As outlined by Roald and Andersson (2018), a number of different options exist for solving this optimization problem under uncertainty. The first one is a one-shot optimization, where the optimal solution is found while respecting all constraints simultaneously (which is demanding in view of the complexity of the problem). The second one is an iterative solution algorithm. The third option similar to the previously mentioned ensemble approach employs Monte Carlo simulations for deriving uncertainty margins and samples uncertainties in power and demand. Here the resulting power flows are calculated for a large number of sample realizations drawn from a given ensemble, leading to Monte Carlo–based uncertainty margins for the voltage. For a comparison of the efficiency of the various approaches, see Roald and Andersson (2018).
IX. Summary and Outlook
Modern power grids are tasked with incorporating increasing shares of renewable energy sources. This ongoing drastic transformation is crucial to mitigating the climate crisis but poses novel challenges to our understanding of the collective dynamics of power grids. Although the underlying equations have been known for a long time, the treatment of such grids on the systems level remains difficult, mainly due to their nonlinearities, distinct heterogeneities in the dynamics of the nodes, and their interaction topology, as well as the various kinds of perturbations and fluctuations present. Additionally, power grids are embedded in and interconnected to several other complex systems, such as energy markets, consumers, and weather. We discussed here central challenges of decentralized and heterogeneous power grids and ways to analyze, understand, design, and influence them.
The review started with a description of power grids based on the principles of physics and the treatment of the most important node models. A crucial aspect here is that renewable sources are connected via electronic inverters. This necessitates an extension of the classic models and leads to fundamental changes in the dynamics of power grids. To study these models and the variety of analysis techniques, datasets on dynamical parameters and network topologies are necessary. However, data on real-world power grids are scarce and often incomplete. Therefore, we presented in Sec. IV typical classes of synthetic models that mimic the main properties of real grids, such as the classic IEEE test cases, but also more recent developments on generating novel network ensembles, which are often closer to real grids and emphasize the physics perspective. Section V discussed an elementary grid containing only one transmission line. This simple bistable system served to introduce the most important static solutions and voltage limits, the bifurcations of the dynamics, and recent stability concepts such as basin stability and stochastic stability, providing a deeper understanding of stability with respect to large and sustained perturbations.
In Sec. VI, recent achievements and insights for realistic networks were presented. To function as a power grid, the networked dynamics needs to have a stable synchronous (phase-locked) state. However, power grids are typically multistable, exhibiting multiple synchronous yet also asynchronous states that are unsuitable for grid operation. This multistability increases substantially if the ubiquitous losses are included, leading to new types of both synchronized and desynchronized states. We gave an overview of the most successful analytic and probabilistic approaches for managing such a rich variety of dynamics, as well as classic linearization-based methods.
The structural stability of grids was the topic of Sec. VII. In complex grids the outage of a single transmission element may induce cascades leading to a large-scale blackout. Such events are well documented in real grids all over the globe, among them several “monster blackouts.” When a single element fails or is added, the network flows reroute in response. A main grid parameter influencing this rerouting is the network distance. But the splitting of a grid into different communities also has a strong impact. New mathematical descriptions of the rerouting process combining methods from statistical physics, nonlinear dynamics, and graph theory are capable of describing these reroutings in new ways and reliably predicting critical links in a network. It turns out that a local stability analysis is often not sufficient and may even be misleading. The occurrence of cascades following individual failures can also be enhanced by transient effects, which shows that dynamical and structural aspects are deeply interwoven. Furthermore, reroutings may yield counterintuitive consequences if transmission elements are added or removed. For instance, under certain conditions the reinforcement of a transition line or the addition of a new line may induce a loss of capacity and stability, a phenomenon known as Braess’s paradox.
The power grid is subject to fluctuations. Section VIII discussed this essential problem with a particular focus on the characteristics introduced by renewables. The non-Gaussian nature of their fluctuations, often characterized by heavy-tailed distributions, such as those due to the cloud structure affecting photovoltaic generation, calls for new modeling and analysis approaches. First, a linear response theory was developed to describe the collective spatiotemporal grid dynamics subject to both stationary and transient and nonstationary and distributed input fluctuations. This led to an efficient way of identifying vulnerability patterns. Second, the emergence of resonances, bulk oscillations, or localizations was explained. Recent research demonstrated using a WKB ansatz that even desynchronization events in grids, which originate from strong but rare fluctuations, can be predicted. The results revise the naive expectation that fat-tail distributions in fluctuating power production always increase the number of rare events of desynchronization. In real networks with renewable components, fluctuations in generation and consumption occur in parallel; i.e., one cannot restrict the robustness analysis to the failure of a single line. Statistical physics provides alternative approaches to treating this problem, such as the consideration of an entire ensemble of power grids, differing in the realization of fluctuations in its microstates, or by chance-constrained ac optimal power flow.
Aspects of power system control and monitoring were only mentioned in this review but will become increasingly important in future power grids with many distributed, fluctuating power sources and low inertia (Milano et al., 2018). This field of research spans from single machines to the entire energy system. Current challenges are the development of control systems for grid-forming inverters that guarantee global dynamic stability (Colombino et al., 2019; Schiffer, Efimov, and Ortega, 2019) or virtual inertia systems that may replace the mechanic inertia of synchronous machines (Chen et al., 2011; Kerdphol et al., 2019). On the grid level, researchers thrive for a better understanding of the interplay of control systems and network dynamics (Tumash, Olmi, and Schöll, 2019; Totz, Olmi, and Schöll, 2020) and the development of new concepts for the control of complex networks (Cornelius, Kath, and Motter, 2013; Huang et al., 2019). Finally, the operation of the large-scale load-frequency control system is subject to multiple external influences, including markets and regulations (Kruse, Schäfer, and Witthaut, 2021a, 2021b).
This review demonstrated that the analysis and design of decentralized power grids demand an interdisciplinary approach where techniques and concepts of power engineering and control theory are newly combined with those from statistical physics, complex dynamical systems, and network science. Such an approach will continue to form the backbone for gaining a deeper understanding of all aspects discussed here and will provide the inspiration for important future research directions.
Future integrated studies will have to connect the power grid to a broader range of important aspects of the energy system, such as the energy market, the information and communication infrastructure driving the control, consumer behavior, electromobility, political and planning constraints, etc. All these aspects point toward the nature of the power grid as part of a multilayered network that is embedded in and connected to other complex systems. It is becoming increasingly clear that machine learning techniques have enormous potential in the study and control of complex heterogeneous systems. Specifically, graph neural networks have vast untapped potential. The question of how to effectively combine such data-driven approaches with physics, engineering-based concepts, and generic modeling remains open in the context of power grids.
Finally, a deeper understanding of fundamental concepts underlying the collective dynamics of power grids should enable more than an energy transition in highly developed countries with well-connected grids required to stop the worsening of the climate crisis. It may also enable an implementation of sustainable power systems in rural areas of the global South and in megacities that now may be designed from scratch. Finally, it should improve an understanding of small-scale energy islands that can play a crucial role in the sustainable electrification of local communities that still remain off grid today.
List of Symbols and Abbreviations
CC OPF chance-constrained optimal power flow DC OPF linearized or dc optimal power flow EMF electromotive force ENTSO-E European Network of Transmission System Operators for Electricity HVDC high voltage directed current LODF line outage distribution factor OPF optimal power flow PI law proportional-integral control law PTDF power transfer distribution factor pu (system) per unit (system of rescaled units) RoCoF rate of change of frequency WKB Wentzel-Kramers-Brillouin
Acknowledgments
This review was produced during our work on the collaborative research project “Kollektive Nichtlineare Dynamik Komplexer Stromnetze” (CoNDyNet), which was funded by the German Federal Ministry of Education and Research (BMBF Grants No. 03EK3055 and No. 03SF0472). D. W. acknowledges support from the Helmholtz Association via the joint initiative “Energy System 2050—A Contribution of the Research Field Energy” and Grant No. VH-NG-1025. J. K. was supported by Russian Ministry of Science and Education Agreement No. 075-15-2020-808.
References
- Agresti, A. and B. A. Coull, Am. Stat. 52, 119 (1998).
- Ainsworth, N. and S. Grijalva, IEEE Trans. Power Syst. 28, 4310 (2013).
- Albert, R., I. Albert, and G. L. Nakarado, Phys. Rev. E 69, 025103 (2004).
- Albert, R., H. Jeong, and A. Barabási, Nature (London) 406, 378 (2000).
- Anderson, P. W., Phys. Rev. 109, 1492 (1958).
- Anvari, M., F. Hellmann, and X. Zhang, Chaos 30, 063140 (2020).
- Anvari, M., G. Lohmann, M. Wächter, P. Milan, E. Lorenz, D. Heinemann, M. R. R. Tabar, and J. Peinke, New J. Phys. 18, 063027 (2016).
- Anvari, M., M. Wächter, and J. Peinke, Europhys. Lett. 116, 60009 (2016).
- Auer, S., F. Hellmann, M. Krause, and J. Kurths, Chaos 27, 127003 (2017).
- Auer, S., K. Kleis, P. Schultz, J. Kurths, and F. Hellmann, Eur. Phys. J. Special Topics 225, 609 (2016).
- Balestra, C., F. Kaiser, D. Manik, and D. Witthaut, Chaos 29, 123119 (2019).
- Bergen, A. R. and D. J. Hill, IEEE Trans. Power Appar. Syst. PAS-100, 25 (1981).
- Berner, R., S. Yanchuk, and E. Schöll, Phys. Rev. E 103, 042315 (2021).
- Bialek, J. W. and V. Vahidinasab, IEEE Trans. Power Syst. 37, 467 (2022).
- Biggs, N., Bull. Lond. Math. Soc. 29, 641 (1997).
- Birchfield, A., https://electricgrids.engr.tamu.edu/, 2016.
- Birchfield, A. B., T. Xu, K. M. Gegner, K. S. Shetye, and T. J. Overbye, IEEE Trans. Power Syst. 32, 3258 (2017).
- Birchfield, A. B., T. Xu, and T. J. Overbye, IEEE Trans. Power Syst. 33, 6667 (2018).
- Bloess, A., W.-P. Schill, and A. Zerrahn, Appl. Energy 212, 1611 (2018).
- Bloomfield, H., D. J. Brayshaw, L. Shaffrey, P. J. Coker, and H. E. Thornton, Environ. Res. Lett. 13, 054028 (2018).
- Bollobás, B., Modern Graph Theory (1998).
- Bompard, E., L. Luo, and E. Pons, Int. J. Crit. Infrastruct. 11, 15 (2015).
- Bosetti, H. and S. Khan, IEEE Trans. Power Syst. 33, 2078 (2018).
- Böttcher, P. C., A. Otto, S. Kettemann, and C. Agert, Chaos 30, 013122 (2020).
- Braess, D., Unternehmensforschung 12, 258 (1968).
- Brummitt, C. D., P. D. Hines, I. Dobson, C. Moore, and R. M. D’Souza, Proc. Natl. Acad. Sci. U.S.A. 110, 12159 (2013).
- Carrasco, J. M., L. G. Franquelo, J. T. Bialasiewicz, E. Galván, R. C. PortilloGuisado, M. M. Prats, J. I. León, and N. Moreno-Alfonso, IEEE Trans. Ind. Electron. 53, 1002 (2006).
- Carreras, B. A., V. E. Lynch, I. Dobson, and D. E. Newman, Chaos 12, 985 (2002).
- Carvalho, R., L. Buzna, R. Gibbens, and F. Kelly, New J. Phys. 17, 095001 (2015).
- Chen, W., D. Wang, J. Liu, T. Basar, K. H. Johansson, and L. Qiu, IFAC-PapersOnLine 49, 97 (2016).
- Chen, Y., R. Hesse, D. Turschner, and H.-P. Beck, in Proceedings of the 2011 International Conference on Power Engineering, Energy and Electrical Drives, Torremolinos, Spain, 2011, (IEEE, 2011), pp. 1–6.
- Chopra, N. and M. W. Spong, IEEE Trans. Autom. Control 54, 353 (2009).
- Chow, Joe H., editor, Power System Coherency and Model Reduction (Springer, New York, 2013).
- Christie, R. D., http://www.ee.washington.edu/research/pstca/, 1999.
- Clauset, A., C. R. Shalizi, and M. E. Newman, SIAM Rev. 51, 661 (2009).
- Coffrin, C., D. Gordon, and P. Scott, 2014 arXiv:1411.0359.
- Coffrin, C., P. Van Hentenryck, and R. Bent, in Proceedings of the 2012 IEEE Power and Energy Society General Meeting, San Diego, 2012, (IEEE, 2012), pp. 1–8.
- Coletta, T., R. Delabays, I. Adagideli, and P. Jacquod, New J. Phys. 18, 103042 (2016).
- Coletta, T. and P. Jacquod, Phys. Rev. E 93, 032222 (2016).
- Coletta, T. and P. Jacquod, IEEE Trans. Control Netw. Syst. 7, 221 (2020).
- Collins, S., P. Deane, B. Ó. Gallachóir, S. Pfenninger, and I. Staffell, Joule 2, 2076 (2018).
- Colombino, M., D. Groß, J.-S. Brouillon, and F. Dörfler, IEEE Trans. Autom. Control 64, 4496 (2019).
- Cornelius, S. P., W. L. Kath, and A. E. Motter, Nat. Commun. 4, 1942 (2013).
- Cotilla-Sanchez, E., P. D. Hines, C. Barrows, and S. Blumsack, IEEE Syst. J. 6, 616 (2012).
- Creutzig, F., P. Agoston, J. Goldschmidt, G. Luderer, G. Nemet, and R. Pietzcker, Nat. Energy 2, 17140 (2017).
- Cuffe, P. and A. Keane, IEEE Syst. J. 11, 1810 (2017).
- Delabays, R., T. Coletta, and P. Jacquod, J. Math. Phys. (N.Y.) 57, 032701 (2016).
- Delabays, R., T. Coletta, and P. Jacquod, J. Math. Phys. (N.Y.) 58, 032703 (2017).
- Dobson, I., B. A. Carreras, V. E. Lynch, and D. E. Newman, Chaos 17, 026103 (2007).
- Dobson, I., B. A. Carreras, and D. E. Newman, Probab. Eng. Inf. Sci. 19, 15 (2005).
- Dörfler, F. and F. Bullo, in Proceedings of the 51st IEEE Conference on Decision and Control (CDC), Maui, 2012, (IEEE, 2012), pp. 7157–7170.
- Dörfler, F. and F. Bullo, SIAM J. Control Optim. 50, 1616 (2012). b.
- Dörfler, F. and F. Bullo, IEEE Trans. Circuits Syst. I 60, 150 (2013). a.
- Dörfler, F. and F. Bullo, in Proceedings of the 2013 IEEE Power & Energy Society General Meeting, Vancouver, British Columbia, Canada, (IEEE, 2013), pp. 1–5.
- Dörfler, F. and F. Bullo, Automatica 50, 1539 (2014).
- Dörfler, F., M. Chertkov, and F. Bullo, Proc. Natl. Acad. Sci. U.S.A. 110, 2005 (2013).
- Dörfler, F., J. W. Simpson-Porco, and F. Bullo, Proc. IEEE 106, 977 (2018).
- Douglass, D. A. and A.-A. Edris, IEEE Trans. Power Delivery 11, 1407 (1996).
- Dykman, M. I., E. Mori, J. Ross, and P. Hunt, J. Chem. Phys. 100, 5735 (1994).
- Edenhofer, O., R. Pichs-Madruga, and Y. Sokona, editors, IPCC Special Report on Renewable Energy Sources and Climate Change Mitigation (Cambridge University Press, Cambridge, England, 2011).
- Engelmann, A., T. Mühlpfordt, Y. Jiang, B. Houska, and T. Faulwasser, IFAC-PapersOnLine 50, 5536 (2017).
- Erseghe, T., IEEE Trans. Power Syst. 29, 2370 (2014).
- Espejo, R., S. Lumbreras, and A. Ramos, IEEE Syst. J. 13, 3050 (2019).
- European Network of Transmission System Operators for Electricity, https://www.entsoe.eu/publications/system-operations-reports/operation-handbook, 2004.
- Fairley, P., IEEE Spectrum 41, 22 (2004).
- Farid, Amro M. et al., https://amfarid.scripts.mit.edu/Datasets/index.html, 2019.
- Fiaz, S., D. Zonetti, R. Ortega, J. M. Scherpen, and A. Van der Schaft, Eur. J. Control 19, 477 (2013).
- Fiedler, M., Czech. Math. J. 23, 298 (1973).
- Figueres, C., H. J. Schellnhuber, G. Whiteman, J. Rockström, A. Hobley, and S. Rahmstorf, Nature (London) 546, 593 (2017).
- Filatrella, G., A. H. Nielsen, and N. F. Pedersen, Eur. Phys. J. B 61, 485 (2008).
- Fleer, J. and P. Stenzel, J. Energy Storage 8, 320 (2016).
- Fliscounakis, S., P. Panciatici, F. Capitanescu, and L. Wehenkel, IEEE Trans. Power Syst. 28, 4909 (2013).
- Frank, S., I. Steponavice, and S. Rebennack, Energy Syst. 3, 221 (2012). a.
- Frank, S., I. Steponavice, and S. Rebennack, Energy Syst. 3, 259 (2012). b.
- Gajduk, A., M. Todorovski, and L. Kocarev, Eur. Phys. J. Special Topics 223, 2387 (2014).
- Gajduk, A., M. Todorovski, J. Kurths, and L. Kocarev, New J. Phys. 16, 115011 (2014).
- Galindo-González, C. C., D. Angulo-García, and G. Osorio, New J. Phys. 22, 103033 (2020).
- Gambuzza, L. V., A. Buscarino, L. Fortuna, M. Porfiri, and M. Frasca, IEEE J. Emerg. Sel. Top. Circuits Syst. 7, 413 (2017).
- García-Mata, I., O. Giraud, B. Georgeot, J. Martin, R. Dubertrand, and G. Lemarié, Phys. Rev. Lett. 118, 166801 (2017).
- Gelbrecht, M., J. Kurths, and F. Hellmann, New J. Phys. 22, 033032 (2020).
- Gorjão, L. R., M. Anvari, H. Kantz, C. Beck, D. Witthaut, M. Timme, and B. Schäfer, IEEE Access 8, 43082 (2020).
- Gorjão, L. R., R. Jumar, H. Maass, V. Hagenmeyer, G. C. Yalcin, J. Kruse, M. Timme, C. Beck, D. Witthaut, and B. Schäfer, Nat. Commun. 11, 6362 (2020).
- Grainger, J. J. and W. D. Stevenson, Jr., Power System Analysis (McGraw-Hill, New York, 1994).
- Gröger, O., H. A. Gasteiger, and J.-P. Suchsland, J. Electrochem. Soc. 162, A2605 (2015).
- Güler, T., G. Gross, and M. Liu, IEEE Trans. Power Syst. 22, 879 (2007).
- Guo, J., Y. Fu, Z. Li, and M. Shahidehpour, IEEE Trans. Power Syst. 24, 1633 (2009).
- Guo, L., C. Liang, A. Zocca, S. H. Low, and A. Wierman, IEEE Trans. Power Syst. 36, 4140 (2021).
- Hackl, J. and B. T. Adey, J. Complex Netw. 7, 254 (2019).
- Haehne, H., K. Schmietendorf, S. Tamrakar, J. Peinke, and S. Kettemann, Phys. Rev. E 99, 050301 (2019).
- Haehne, H., J. Schottler, M. Waechter, J. Peinke, and O. Kamps, Europhys. Lett. 121, 30001 (2018).
- Han, C., D. Witthaut, M. Timme, and M. Schröder, PLoS One 14, e0225346 (2019).
- Hata, S. and H. Nakao, Sci. Rep. 7, 1121 (2017).
- Heide, D., L. von Bremen, M. Greiner, C. Hoffmann, M. Speckmann, and S. Bofinger, Renewable Energy 35, 2483 (2010).
- Hellmann, F., P. Schultz, C. Grabow, J. Heitzig, and J. Kurths, Sci. Rep. 6, 29654 (2016).
- Hellmann, F., P. Schultz, P. Jaros, R. Levchenko, T. Kapitaniak, J. Kurths, and Y. Maistrenko, Nat. Commun. 11, 592 (2020).
- Hindes, J., P. Jacquod, and I. B. Schwartz, Phys. Rev. E 100, 052314 (2019).
- Hines, P., E. Cotilla-Sanchez, and S. Blumsack, Chaos 20, 033122 (2010).
- Hines, P. D., I. Dobson, and P. Rezaei, IEEE Trans. Power Syst. 32, 958 (2017).
- Hörsch, J., F. Hofmann, D. Schlachtberger, and T. Brown, Energy Strategy Rev. 22, 207 (2018).
- Hörsch, J., H. Ronellenfitsch, D. Witthaut, and T. Brown, Electr. Power Syst. Res. 158, 126 (2018).
- Huang, L., J. Coulson, J. Lygeros, and F. Dörfler, in Proceedings of the 58th IEEE Conference on Decision and Control (CDC), (IEEE, 2019), pp. 8130–8135.
- Huneault, M. and F. D. Galiana, IEEE Trans. Power Syst. 6, 762 (1991).
- Hutcheon, N. and J. Bialek, in Proceedings of the IEEE PowerTech Conference, Grenoble, 2013, (IEEE, 2013).
- IPCC, Climate Change 2014: Mitigation of Climate Change. Contribution of Working Group 3 to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (Cambridge University Press, Cambridge, England, 2014).
- Jadbabaie, A., N. Motee, and M. Barahona, in Proceedings of the 2004 American Control Conference, Boston, 2004, No. 5 (IEEE, 2004), pp. 4296–4301.
- Jafarpour, S., E. Y. Huang, K. D. Smith, and F. Bullo, 2019 arXiv:1901.11189.
- Jung, D. and S. Kettemann, Phys. Rev. E 94, 012307 (2016).
- Kaiser, F., V. Latora, and D. Witthaut, Nat. Commun. 12, 3143 (2021).
- Kaiser, F., H. Ronellenfitsch, V. Latora, and D. Witthaut, 2021 arXiv:2105.06687.
- Kaiser, F., H. Ronellenfitsch, and D. Witthaut, Nat. Commun. 11, 5796 (2020).
- Kaiser, F., J. Strake, and D. Witthaut, New J. Phys. 22, 013053 (2020).
- Kaiser, F. and D. Witthaut, Phys. Rev. Research 3, 023161 (2021). a.
- Kaiser, F. and D. Witthaut, IEEE Access 9, 67364 (2021). b.
- Kashima, K., H. Aoyama, and Y. Ohta, in Proceedings of the 54th IEEE Conference on Decision and Control (CDC), Osaka, Japan, 2015, (IEEE, 2015), pp. 1852–1857.
- Katifori, E., G. J. Szöllősi, and M. O. Magnasco, Phys. Rev. Lett. 104, 048704 (2010).
- Kavasseri, R. and C. Ababei, http://www.dejazzer.com/reds.html, 2021.
- Kerdphol, T., F. S. Rahman, M. Watanabe, Y. Mitani, D. Turschner, and H.-P. Beck, IEEE Access 7, 14422 (2019).
- Kettemann, S., Phys. Rev. E 94, 062311 (2016).
- Kim, H., M. J. Lee, S. H. Lee, and S.-W. Son, Chaos 29, 103132 (2019).
- Kim, H., S. H. Lee, J. Davidsen, and S.-W. Son, New J. Phys. 20, 113006 (2018).
- Kim, H., S. H. Lee, and P. Holme, New J. Phys. 17, 113005 (2015).
- Kim, H., S. H. Lee, and P. Holme, Phys. Rev. E 93, 062318 (2016).
- Klein, D. J. and M. Randić, J. Math. Chem. 12, 81 (1993).
- Klein, M., G. J. Rogers, and P. Kundur, IEEE Trans. Power Syst. 6, 914 (1991).
- Korjani, S., A. Facchini, M. Mureddu, G. Caldarelli, and A. Damiano, Sci. Rep. 8, 16658 (2018).
- Korkali, M., J. G. Veneman, B. F. Tivnan, J. P. Bagrow, and P. D. Hines, Sci. Rep. 7, 44499 (2017).
- Kruse, J., B. Schäfer, and D. Witthaut, Patterns 2, 100365 (2021). a.
- Kruse, J., B. Schäfer, and D. Witthaut, 2021b, arXiv:2109.04802.
- Kuehn, C. and S. Throm, SIAM J. Appl. Math. 79, 1271 (2019).
- Kundur, P., Power System Stability and Control (McGraw-Hill, New York, 1994).
- Kuramoto, Y., in International Symposium on on Mathematical Problems in Theoretical Physics edited by H. Araki (Springer, New York, 1975), p. 420.
- Labavic, D., R. Suciu, H. Meyer-Ortmanns, and S. Kettemann, Eur. Phys. J. Special Topics 223, 2517 (2014).
- Liang, X. and C. Andalib-Bin-Karim, IEEE Trans. Ind. Appl. 54, 3100 (2018).
- Lindner, M. and F. Hellmann, Phys. Rev. E 100, 022124 (2019).
- Liu, H., Z. Hu, Y. Song, and J. Lin, IEEE Trans. Power Syst. 28, 3480 (2013).
- Liu, Z., X. He, Z. Ding, and Z. Zhang, IEEE Trans. Ind. Inf. 15, 1450 (2019).
- Lu, Ning, Ruisheng Diao, R. P. Hafen, N. Samaan, and Y. V. Makarov, in Proceedings of the 2013 IEEE Power & Energy Society General Meeting, Vancouver, British Columbia, Canada, 2013, (IEEE, 2013), pp. 1–5.
- Lun, I. Y. and J. C. Lam, Renewable Energy 20, 145 (2000).
- Ma, J., Y. Sun, X. Yuan, J. Kurths, and M. Zhan, PLoS One 11, e0165943 (2016).
- Machowski, J., J. Bialek, and J. Bumby, Power System Dynamics, Stability and Control (John Wiley & Sons, New York, 2008).
- Manik, D., M. Rohden, H. Ronellenfitsch, X. Zhang, S. Hallerberg, D. Witthaut, and M. Timme, Phys. Rev. E 95, 012319 (2017).
- Manik, D., M. Timme, and D. Witthaut, Chaos 27, 083123 (2017).
- Manik, D., D. Witthaut, B. Schäfer, M. Matthiae, A. Sorge, M. Rohden, E. Katifori, and M. Timme, Eur. Phys. J. Special Topics 223, 2527 (2014).
- Medjroubi, W., C. Matke, and D. Kleinhans, https://www.scigrid.de/, 2015.
- Mehrmann, V., R. Morandin, S. Olmi, and E. Schöll, Chaos 28, 101102 (2018).
- Mehta, D., N. S. Daleo, F. Dörfler, and J. D. Hauenstein, Chaos 25, 053103 (2015).
- Mehta, D., D. K. Molzahn, and K. Turitsyn, in Proceedings of the 2016 American Control Conference (ACC), Boston, 2016, (IEEE, 2016), pp. 1753–1765.
- Menck, P. J., J. Heitzig, J. Kurths, and H. J. Schellnhuber, Nat. Commun. 5, 3969 (2014).
- Menck, P. J., J. Heitzig, N. Marwan, and J. Kurths, Nat. Phys. 9, 89 (2013).
- Milan, P., M. Wächter, and J. Peinke, Phys. Rev. Lett. 110, 138701 (2013).
- Milano, F., Power System Modelling and Scripting (Springer Science+Business Media, New York, 2010).
- Milano, F., F. Dörfler, G. Hug, D. J. Hill, and G. Verbič, in Proceedings of the 2018 Power Systems Computation Conference (PSCC), Dublin, 2018, (IEEE, 2018), pp. 1–25.
- Momoh, J. A., R. Adapa, and M. El-Hawary, IEEE Trans. Power Syst. 14, 96 (1999).
- Momoh, J. A., M. El-Hawary, and R. Adapa, IEEE Trans. Power Syst. 14, 105 (1999).
- Motter, A. E., Phys. Rev. Lett. 93, 098701 (2004).
- Motter, A. E. and Y.-C. Lai, Phys. Rev. E 66, 065102 (2002).
- Motter, A. E., S. A. Myers, M. Anghel, and T. Nishikawa, Nat. Phys. 9, 191 (2013).
- Mureddu, M., G. Caldarelli, A. Chessa, A. Scala, and A. Damiano, PLoS One 10, e0135312 (2015).
- Mureddu, M. and A. Damiano, in Proceedings of the 26th IEEE International Symposium on Industrial Electronics (ISIE), Edinburgh, 2017, (IEEE, 2017), pp. 2069–2074.
- Mureddu, M. and H. Meyer-Ortmanns, Physica (Amsterdam) 490A, 1324 (2018).
- Nesti, T., F. Sloothaak, and B. Zwart, Phys. Rev. Lett. 125, 058301 (2020).
- Nesti, T., A. Zocca, and B. Zwart, Phys. Rev. Lett. 120, 258301 (2018).
- Newman, M. E., Nat. Phys. 8, 25 (2012).
- Newman, M. E. J., Networks—An Introduction (Oxford University Press, Oxford, 2010).
- Nicolaisen, J., V. Petrov, and L. Tesfatsion, IEEE Trans. Evol. Comput. 5, 504 (2001).
- Nishikawa, T. and A. E. Motter, Phys. Rev. E 73, 065106 (2006).
- Nishikawa, T. and A. E. Motter, New J. Phys. 17, 015012 (2015).
- Nitzbon, J., P. Schultz, J. Heitzig, J. Kurths, and F. Hellmann, New J. Phys. 19, 033029 (2017).
- Ochab, J. and P. Gora, Acta Phys. Pol. B Proc. Suppl. 3, 453 (2010).
- Ódor, G. and B. Hartmann, Phys. Rev. E 98, 022305 (2018).
- Oeding, D. and B. R. Oswald, Elektrische Kraftwerke Und Netze (Springer, New York, 2016).
- Olmi, S., Chaos 25, 123125 (2015).
- Olmi, S., A. Navas, S. Boccaletti, and A. Torcini, Phys. Rev. E 90, 042905 (2014).
- Open Energy Modelling Initiative, http://www.openmod-initiative.org/, 2017.
- Open Power Systems Data, http://open-power-system-data.org/, 2017.
- Pacific Northwest National Laboratory and National Rural Electric Cooperative Association, and https://egriddata.org/, 2017.
- Pagnier, L. and P. Jacquod, PLoS One 14, e0213550 (2019). a.
- Pagnier, L. and P. Jacquod, IEEE Access 7, 145889 (2019). b.
- Pagnier, L. and P. Jacquod, c, 10.5281/zenodo.2642175, 2019 http://dx.doi.org/10.5281/zenodo.2642175.
- Parks, P. C., IMA J. Math. Control Inf. 9, 275 (1992).
- Pesch, T., H.-J. Allelein, and J.-F. Hake, Eur. Phys. J. Special Topics 223, 2561 (2014).
- Plietzsch, A., S. Auer, J. Kurths, and F. Hellmann, 2019 arXiv:1903.09585.
- Plietzsch, A., P. Schultz, J. Heitzig, and J. Kurths, Eur. Phys. J. Special Topics 225, 551 (2016).
- Pourbeik, P., P. Kundur, and C. Taylor, IEEE Power Energy Mag. 4, 22 (2006).
- Purchala, K., L. Meeus, D. Van Dommelen, and R. Belmans, in Proceedings of the IEEE Power Engineering Society General Meeting, San Francisco, 2005, (IEEE, 2005).
- Reseau de Transport d’Electricite, https://clients.rte-france.com/lang/an/visiteurs/vie/vie_frequence.jsp (accessed February 2018), 2017.
- Risken, H., The Fokker-Planck Equation (Springer, Berlin, 1996).
- Ritmeester, T. and H. Meyer-Ortmanns, Physica (Amsterdam) 573A, 125927 (2021).
- Roald, L. and G. Andersson, IEEE Trans. Power Syst. 33, 2906 (2018).
- Rockström, J., O. Gaffney, J. Rogelj, M. Meinshausen, N. Nakicenovic, and H. J. Schellnhuber, Science 355, 1269 (2017).
- Rodrigues, F. A., T. K. D. Peron, P. Ji, and J. Kurths, Phys. Rep. 610, 1 (2016).
- Rogelj, J., G. Luderer, R. Pietzcker, E. Kriegler, M. Schaeffer, M. Krey, and K. Riahi, Nat. Clim. Change 5, 519 (2015).
- Rogers, G., Power System Oscillations (Springer Science+Business Media, New York, 2012).
- Rohden, M., D. Jung, S. Tamrakar, and S. Kettemann, Phys. Rev. E 94, 032209 (2016).
- Rohden, M., A. Sorge, M. Timme, and D. Witthaut, Phys. Rev. Lett. 109, 064101 (2012).
- Rohden, M., D. Witthaut, M. Timme, and H. Meyer-Ortmanns, New J. Phys. 19, 013002 (2017).
- Ronellenfitsch, H., D. Manik, J. Horsch, T. Brown, and D. Witthaut, IEEE Trans. Power Syst. 32, 4060 (2017).
- Sauer, P. W., M. A. Pai, and J. H. Chow, Power System Dynamics and Stability: With Synchrophasor Measurement and Power System Toolbox (John Wiley & Sons, Hoboken, NJ, 2018).
- Schäfer, B., C. Beck, K. Aihara, D. Witthaut, and M. Timme, Nat. Energy 3, 119 (2018).
- Schäfer, B., M. Matthiae, X. Zhang, M. Rohden, M. Timme, and D. Witthaut, Phys. Rev. E 95, 060203 (2017).
- Schäfer, B., M. Timme, and D. Witthaut, in Proceedings of the 2018 IEEE PES Innovative Smart Grid Technologies Conference Europe (ISGT-Europe), Sarajevo, 2018, (IEEE, 2018), pp. 1–5.
- Schäfer, B., D. Witthaut, M. Timme, and V. Latora, Nat. Commun. 9, 1975 (2018).
- Schiffer, J., D. Efimov, and R. Ortega, in Proceedings of the 57th IEEE Conference on Decision and Control (CDC), Miami, 2018, (IEEE, 2018), pp. 800–805.
- Schiffer, J., D. Efimov, and R. Ortega, Automatica 109, 108550 (2019).
- Schiffer, J., R. Ortega, A. Astolfi, J. Raisch, and T. Sezi, Automatica 50, 2457 (2014).
- Schlott, M., A. Kies, T. Brown, S. Schramm, and M. Greiner, Appl. Energy 230, 1645 (2018).
- Schmietendorf, K., J. Peinke, R. Friedrich, and O. Kamps, Eur. Phys. J. Special Topics 223, 2577 (2014).
- Schmietendorf, K., J. Peinke, and O. Kamps, Eur. Phys. J. B 90, 222 (2017).
- Schneider, K. P. et al., http://sites.ieee.org/pes-testfeeders/, 1991.
- Schneider, K. P. et al., IEEE Trans. Power Syst. 33, 3181 (2018).
- Schröder, M., M. Timme, and D. Witthaut, Chaos 27, 073119 (2017).
- Schultz, P., J. Heitzig, and J. Kurths, New J. Phys. 16, 125001 (2014). a.
- Schultz, P., J. Heitzig, and J. Kurths, Eur. Phys. J. Special Topics 223, 2593 (2014). b.
- Schultz, P., F. Hellmann, J. Heitzig, and J. Kurths, 2016 arXiv:1701.06968.
- Schultz, P., F. Hellmann, K. N. Webster, and J. Kurths, Chaos 28, 043102 (2018).
- Seguro, J. and T. Lambert, J. Wind Eng. Ind. Aerodyn. 85, 75 (2000).
- Semerow, A., S. Höhn, M. Luther, W. Sattinger, H. Abildgaard, A. D. Garcia, and G. Giannuzzi, in Proceedings of the 2015 IEEE PowerTech Conference, (IEEE, 2015).
- Setreus, J. and L. Bertling, in Proceedings of the 10th International Conference on Probablistic Methods Applied to Power Systems, Rincon, PR, 2008, (IEEE, 2008), pp. 1–8.
- Sharafutdinov, K., L. Rydin Gorjão, M. Matthiae, T. Faulwasser, and D. Witthaut, Chaos 28, 033117 (2018).
- Simonsen, I., L. Buzna, K. Peters, S. Bornholdt, and D. Helbing, Phys. Rev. Lett. 100, 218701 (2008).
- Simpson-Porco, J. W., IEEE Trans. Power Syst. 33, 2477 (2018).
- Skar, S. J., SIAM J. Appl. Math. 39, 475 (1980).
- Solé, R. V., M. Rosas-Casals, B. Corominas-Murtra, and S. Valverde, Phys. Rev. E 77, 026102 (2008).
- Soltan, S., A. Loh, and G. Zussman, IEEE Syst. J. 13, 625 (2019).
- Soltan, S. and G. Zussman, in Proceedings of the 2016 IEEE Power and Energy Society General Meeting (PESGM), Boston, 2016, (IEEE, 2016), pp. 1–5.
- Song, Y., D. J. Hill, and T. Liu, in Proceedings of the 2015 IEEE Conference on Control Applications (CCA), Sydney, 2015, (IEEE, 2015), p. 201.
- Staffell, I. and S. Pfenninger, Energy 145, 65 (2018).
- Stott, B., Proc. IEEE 67, 219 (1979).
- Stott, B., J. Jardim, and O. Alsac, IEEE Trans. Power Syst. 24, 1290 (2009).
- Strake, J., F. Kaiser, F. Basiri, H. Ronellenfitsch, and D. Witthaut, New J. Phys. 21, 053009 (2019).
- Taher, H., S. Olmi, and E. Schöll, Phys. Rev. E 100, 062306 (2019).
- Tamrakar, S., M. Conrath, and S. Kettemann, Sci. Rep. 8, 6459 (2018).
- Taylor, R., J. Phys. A 45, 055102 (2012).
- Torres-Sánchez, L., G. T. Freitas de Abreu, and S. Kettemann, Phys. Rev. E 101, 012313 (2020).
- Totz, C. H., S. Olmi, and E. Schöll, Phys. Rev. E 102, 022311 (2020).
- Trias, A., in Proceedings of the 2012 Power and Energy Society General Meeting, San Diego, 2012, (IEEE, 2012), pp. 1–8.
- Tumash, L., S. Olmi, and E. Schöll, Europhys. Lett. 123, 20001 (2018).
- Tumash, L., S. Olmi, and E. Schöll, Chaos 29, 123105 (2019).
- Tyloo, M., T. Coletta, and P. Jacquod, Phys. Rev. Lett. 120, 084101 (2018).
- Tyloo, M. and P. Jacquod, IEEE Control Syst. Lett. 5, 929 (2021).
- Tyloo, M., L. Pagnier, and P. Jacquod, Sci. Adv. 5, eaaw8359 (2019).
- UCTE Operations Handbook, https://www.entsoe.eu/fileadmin/user_upload/_library/publications/entsoe/Operation_Handbook/introduction_v25.pdf, 2009.
- Union for the Co-ordination of Transmission of Electricity, https://www.entsoe.eu/fileadmin/user_upload/_library/publications/ce/otherreports/20040427_UCTE_IC_Final_report.pdf, 2004.
- Union for the Coordination of Transmission of Electricity, https://www.entsoe.eu/fileadmin/user_upload/_library/publications/ce/otherreports/Final-Report-20070130.pdf, 2007.
- U.S.-Canada Power System Outage Task Force, https://energy.gov/sites/prod/files/oeprod/DocumentsandMedia/BlackoutFinal-Web.pdf, 2014.
- Van Cutsem, T. and C. Vournas, Voltage Stability of Electric Power Systems (Springer Science+Business Media, New York, 2007).
- Van der Vleuten, E. and V. Lagendijk, Energy Policy 38, 2042 (2010).
- Van Hertem, D., J. Verboomen, K. Purchala, R. Belmans, and W. Kling, in Proceedings of the 8th IEE International Conference on AC and DC Power Transmission, London, 2006, (IEEE, 2006), pp. 58–62.
- van Kampen, N., Stochastic Processes in Physics and Chemistry (Elsevier, Amsterdam, 2007).
- Venkatasubramanian, V. and Y. Li, in Proceedings of Bulk Power System Dynamics and Control–VI, Cortina d’Ampezzo, Italy, 2004, (IEEE, 2004), p. 685.
- Verboomen, J., D. Van Hertem, P. H. Schavemaker, W. L. Kling, and R. Belmans, in Proceedings of the 2005 International Conference on Future Power Systems, Amsterdam, 2005, (IEEE, 2005), p. 6.
- Villella, F., S. Leclerc, I. Erlich, and S. Rapoport, in Proceedings of the 3rd IEEE PES Innovative Smart Grid Technologies Europe (ISGT Europe), Berlin, 2012, (IEEE, 2012), pp. 1–9.
- Wang, Z., A. Scaglione, and R. J. Thomas, IEEE Trans. Smart Grid 1, 28 (2010).
- Weber, J., J. Wohland, M. Reyers, J. Moemken, C. Hoppe, J. G. Pinto, and D. Witthaut, PLoS One 13, e0201457 (2018).
- Weckesser, T., H. Jóhannsson, and J. Østergaard, in Proceedings of the 2013 IREP Symposium Bulk Power System Dynamics and Control–IX, Rethymno, Greece, 2013, (IEEE, 2013), pp. 1–9.
- Wiedermann, M., J. F. Donges, J. Kurths, and R. V. Donner, Phys. Rev. E 93, 042308 (2016).
- Wiegmans, B., 10.5281/zenodo.55853, 2016 http://dx.doi.org/10.5281/zenodo.55853.
- Wiser, R., K. Jenni, J. Seel, E. Baker, M. Hand, E. Lantz, and A. Smith, Nat. Energy 1, 16135 (2016).
- Witthaut, D., M. Rohden, X. Zhang, S. Hallerberg, and M. Timme, Phys. Rev. Lett. 116, 138701 (2016).
- Witthaut, D. and M. Timme, New J. Phys. 14, 083036 (2012).
- Witthaut, D. and M. Timme, Phys. Rev. E 92, 032809 (2015).
- Wohland, J., N. E. Omrani, N. Keenlyside, and D. Witthaut, Wind Energy Sci. 4, 515 (2019).
- Wohland, J., M. Reyers, J. Weber, and D. Witthaut, Earth Syst. Dyn. 8, 1047 (2017).
- Wolff, M. F., P. G. Lind, and P. Maass, Chaos 28, 103120 (2018).
- Wolff, M. F., K. Schmietendorf, P. G. Lind, O. Kamps, J. Peinke, and P. Maass, Chaos 29, 103149 (2019).
- Wood, A. J., B. F. Wollenberg, and G. B. Sheblé, Power Generation, Operation and Control (John Wiley & Sons, New York, 2014).
- Xi, K., J. L. Dubbeldam, and H. X. Lin, Chaos 27, 013109 (2017).
- Xu, T., A. B. Birchfield, and T. J. Overbye, IEEE Trans. Power Syst. 33, 6501 (2018).
- Yang, Y. and A. E. Motter, Phys. Rev. Lett. 119, 248302 (2017).
- Yang, Y., T. Nishikawa, and A. E. Motter, Science 358, eaan3184 (2017).
- Yang, Z., K. Xie, J. Yu, H. Zhong, N. Zhang, and Q. Xia, IEEE Trans. Power Syst. 34, 1315 (2019).
- You, H., V. Vittal, and X. Wang, IEEE Trans. Power Syst. 19, 483 (2004).
- Zhang, H., G. T. Heydt, V. Vittal, and J. Quintero, IEEE Trans. Power Syst. 28, 3471 (2013).
- Zhang, H. and P. Li, IET Gener. Transm. Distrib. 4, 553 (2010).
- Zhang, X., S. Hallerberg, M. Matthiae, D. Witthaut, and M. Timme, Sci. Adv. 5, eaav1027 (2019).
- Zhang, X., D. Witthaut, and M. Timme, Phys. Rev. Lett. 125, 218301 (2020).
- Zhang, X.-P., C. Rehtanz, and B. Pal, Flexible AC Transmission Systems: Modelling and Control (Springer, Heidelberg, 2012).
- Zhang, Y., L. Wehenkel, P. Rousseaux, and M. Pavella, Int. J. Electr. Power Energy Syst. 19, 195 (1997).
- Zhou, K., I. Dobson, Z. Wang, A. Roitershtein, and A. P. Ghosh, IEEE Trans. Power Syst. 35, 3224 (2020).
- Zimmerman, R. D., C. E. Murillo-Sanchez, and R. J. Thomas, IEEE Trans. Power Syst. 26, 12 (2011).
- Zocca, A., C. Liang, L. Guo, S. H. Low, and A. Wierman, 2021 arXiv:2105.05234.
Issue
Operations in the APS Offices, including the Editorial Office, will pause starting Friday afternoon, December 23, 2022 through Monday, January 2, 2023. Journal articles will continue to be published December 23 - 30, 2022. No articles will be published on January 2, 2023. Submissions, referee reports, and other correspondence will be received and timestamped for processing. Normal business operations will resume on Tuesday, January 3, 2023. We appreciate your understanding as processing and response times will be delayed.

